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1.  Introduction 

Cruise  missiles  over  land  and  sea  cluttered  background  are  serious  threats  to  Infrared  Search 
and  Track  Systems  (IRST’s).  In  general,  these  threats  are  stealth  in  both  the  infrared  and  radio 
frequency  bands.  That  is,  their  thermal  infrared  signature  and  their  radar  cross  section  can  be 
quite  small.  Future  predicted  threats,  i.e.  the  next  generation  of  cruise  missiles,  will  be  even 
more  difficult  to  detect  at  a  sufficient  range  to  counter.  Further,  low  elevation  trajectory  objects, 
such  as  sea  skimming  missiles,  have  radar  signals  with  large  amounts  of  temporally  and  spatially 
correlated  interference  called  multipath.  This  multipath  problem  remains  an  enormous  obstacle 
to  existing  trackers.  Hence,  new  technology  is  needed  which  will  allow  for  the  timely  detection, 
tracking,  and  identification  of  such  threats. 

IRST  systems  are  one  component  of  a  multisensor  suite  which  can  meet  the  technical  challenge 
of  the  timely  detection/track/identification  of  low  signal-to-(noise+clutter)  ratio  (S(N+C)R) 
targets.  The  multisensor  suite  should  include  an  IRST,  Radar,  and  a  coherent  laser  (Lidar).  We 
envision  a  cueing  hierarchy  where  the  IRST  can  cue  the  Radar  or  the  Radar  cues  the  IRST.  Once 
a  candidate  track  is  established  the  Lidar  can  be  used  to  identify  the  target  by  its  micro  doppler 
signature. 

In  this  report  we  describe  the  developed  computationally  efficient  algorithms  and  adaptive  archi¬ 
tecture  with  optimized  overall  performance  (statistical  and  computational)  for  real-time  reliable 
detection  and  tracking  of  low-observable  targets  in  IRST  systems.  Despite  the  fact  that  we  fo¬ 
cus  on  an  IRST  against  cruise  missiles  over  land  and  sea  cluttered  backgrounds,  the  results  are 
equally  applicable  to  other  sensors  (e.g.,  Radar,  Lidar). 

In  the  research  we  concentrated  on  the  three  interrelated  problems:  (1)  efficient  clutter  suppres¬ 
sion;  (2)  development  of  the  adaptive  track-before-detect  (TbD)  architecture  based  on  optimal 
nonlinear  filtering  (ONF);  (2)  development  of  efficient  algorithms  for  detection  of  a  priori  un¬ 
known  number  of  targets  that  may  appear  and  disappear  at  unknown  points  in  time. 

\ 

The  report  is  organized  as  follows.  In  Section  2  we  formulate  the  problem,  describe  a  structure  of 
the  system  to  be  developed,  outline  popular  track-before-detect  methods  that  are  in  current  use 
and  suggest  an  alternative  method,  which  is  based  on  the  optimal  nonlinear  filtering.  In  Section  3 
we  describe  basic  models  and  assumptions  on  signals  and  clutter  that  are  used  in  developing  of 
signal  processing  algorithms  (clutter  suppression,  track-before-detect)  and  detection  algorithms. 
In  Section  4,  two  clutter  removal  algorithms  are  presented  based  on  nonparametric  and  semi- 
parametric  spatio-temporal  filtering.  Section  5  is  especially  important.  Here  we  describe  an 
optimal  nonlinear  filtering  technique  that  is  used  for  track-before-detect  of  very  low  observable 
targets.  Also,  we  show  that  the  proposed  method  is  a  complete  generalization  of  the  multi¬ 
dimensional  (spatio-temporal)  matched  filtering  (particularly,  3D  matched  filter).  Furthermore 
it  is  shown  that  the  spatio-temporal  matched  filter  coincides  with  the  developed  optimal  nonlin¬ 
ear  filter  when  a  target  moves  according  to  deterministic  trajectory.  In  Section  6  we  develop  three 
track  appearance/disappearance  algorithms  all  of  which  take  into  account  requirements  specific 
for  surveillance  systems.  The  results  of  simulation  and  processing  real  IR  data  obtained  from 
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SPAWAR  (Space  &  Naval  Warfare  Systems  Center,  San  Diego,  CA)  are  presented  in  Section  7. 
Finally  in  Section  8  we  provide  a  conclusion  and  the  plan  of  our  future  research. 


2.  Problem  Formulation  and  Background 


The  generalized  block-diagram  of  the  system  under  investigation  is  shown  in  Figure  1.  We 
develop  both  the  signal  processing  architecture  (clutter  removal  algorithms  and  TbD  algorithms) 
and  track  detection  algorithms. 


_ SignalProcessing JBlock _ 


Detections  (Blips) 


Figure  1.  Generalized  block-diagram  of  the  developed  system 


New  algorithms  concerning  the  stages  of  data  processing  specified  above  are  developed  under  the 
following  realistic  conditions. 

•  Cluttered  background  is  much  more  intensive  than  both  equivalent  intrinsic  (instrumental) 
noise  of  the  sensor  and  signal  intensity  of  the  targets  to  be  detected.  This  causes  a  necessity 
of  practically  complete  suppression  of  a  clutter. 

•  Exterior  conditions  of  observation  are  characterized  by  an  extremely  high  variability  and 
prior  uncertainty  and  may  not  be  predicted  with  sufficient  accuracy. 

•  Prior  information  that  is  needed  to  develop  ideal  (Bayes)  data  processing  algorithms  is  not 
available.  Particularly,  statistical  models  of  signals  and  moreover  exterior  background  as 
well  as  the  models  of  changing  target  situation  are  extremely  unreliable.  Such  models  can  be 
useful  only  as  tools  for  performance  evaluation  in  certain  scenarios  but  not  for  development 
of  data  processing  algorithms.  Practical  data  processing  algorithms  should  be  developed 
on  the  base  of  sequential  application  of  robust,  adaptive  and  minimax  methods  that  are 
invariant  to  the  prior  uncertainty. 

•  In  practice  the  estimated  parameters  of  targets  (for  example,  trajectory  parameters)  should 
be  guaranteed  for  any  degree  of  prior  uncertainty.  Specifically,  each  estimate  should  be 
supplemented  with  an  appropriate  domain  of  minimum  size  that  contains  the  real  unknown 
value  of  estimated  parameter  with  100%  (or  at  least  (1  —  e)%)  assurance. 
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2.1.  TbD  Methods:  Merits  and  Drawbacks.  The  most  challenging  problem  for  an  IRST 
system  is  the  detection  of  a  maneuvering  target  against  a  strong  clutter  background.  To  illustrate 
the  importance  of  this  task,  we  remark  that  under  certain  conditions  a  few  seconds  decrease  in  the 
time  it  takes  to  detect  a  sea/surface  skimming  cruise  missile  can  yield  a  significant  increase  in  the 
probability  of  raid  annihilation.  The  problem  of  detection  is  extremely  difficult  in  low  S(N+C)R 
when  localization  of  the  target  based  upon  a  single  non-stationary  image  is  impossible.  In  this 
case  one  has  to  align  successive  frames  according  to  typical  patterns  of  target  dynamics  and  any 
results  of  “preliminary”  tracking.  This  approach  to  detection  of  a  low  S(N+C)R  target  is  usually 
referred  to  as  “track-before-detect”  (TbD).  Its  success  depends  crucially  on  the  quality  of  the 
“preliminary”  tracker.  Thus,  the  development  of  the  efficient  coherent  signal  processing  based 
on  the  TbD  methodology  becomes  crucial  point  in  the  low-observable  target  detection  problem. 
In  contrast  to  other  TbD  methods,  we  solve  this  problem  by  applying  optimal  nonlinear  filtering 
approach. 

We  now  overview  several  popular  methods  for  tracking  before  detection  that  are  in  current  use. 


2.1.1.  Spatio-Temporal  Matched  Filters  and  Velocity  Filter  Banks.  Given  a  sequence 
of  frames,  consider  the  problem  of  detecting  a  small  (unresolved)  target  against  a  clutter  back¬ 
ground.  The  probabilities  of  errors  can  be  decreased  by  applying  a  3-D  (spatio-temporal) 
matched  filter  prior  to  detection  (thresholding)  operation.  If  the  spatial  distribution  of  the  tar¬ 
get  and  its  velocity  remain  unchanged  (i.e.,  if  targets  move  with  known  constant  speed  along 
a  line  on  the  plane)  and  the  noise  and  the  clutter  are  Gaussian  processes,  the  3-D  filter  is 
the  optimal  method  of  detection,  see  [38].  It  is  easily  shown  that  the  same  result  is  true  for 
more  dimensions,  i.e.  a  (m  +  1)-D  matched  filter  is  an  optimal  method  of  processing  under 
aforementioned  conditions  ( m  is  the  spatial  dimensionality).  Also,  under  these  conditions  the 
“target”  component  of  the  multi-dimensional  matched  filter  separates  into  spatial  and  temporal 
components.  The  temporal  component  is  usually  called  a  velocity  filter. 

Typically  the  target  velocity  is  unknown  and  hence  the  single  velocity  filter  cannot  be  used.  This 
problem  is  usually  overcome  by  hypothesizing  velocities  and  implementing  a  velocity  filter  bank 
(see,  e.g.,  [43,  48,  50]). 

One  of  the  main  drawbacks  of  spatio-temporal  matched  filters  and  other  modifications  such  as 
banks  of  assumed  velocity  filters  is  that  they  are  not  able  to  work  with  maneuvering  targets. 
Performance  of  the  algorithms  is  substantially  degraded  in  the  presence  of  velocity  mismatch  or 
in  event  of  target  maneuver. 


2.1.2.  Dynamic  Programming  Methods.  The  Dynamic  Programming  methods  for  TbD 
showed  a  big  advantage  over  the  conventional  MHT  method  and  over  the  3— D  matched  fil¬ 
ter  with  velocity  mismatch  [1,  2,  7,  21,  57].  Particularly  the  results  of  Fernandez  et  al.  [21]  show 
that  application  of  the  Viterbi  TbD  algorithm  over  10  frames  of  IR  data  yields  about  a  7  dB 
improvement  in  detection  sensitivity  over  conventional  thresholding/peak-detection  procedures. 
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This  approach  avoided  problems  with  velocity  mismatch  and  could  handle  targets  with  slow 
maneuvers.  However,  good  performance  of  dynamic  programming  algorithms  is  observed  for 
moderate  S(N+C)R  (over  3  dB  after  preprocessing  and  clutter  suppression)  with  rapid  degra¬ 
dation  as  S(N+C)R  reduces  further  [57].  In  addition,  the  computational  complexity  of  these 
sophisticated  methods  is  fairly  high. 


2.1.3.  Extended  Kalman  Filter.  To  date,  the  extended  Kalman  filter  (EKF)  along  with  minor 
variations  have  been  the  dominant  algorithm  technology  in  real-time  tracking.  It  is  the  basis 
for  current  practical  single-  or  multi-target  trackers  for  point  equivalent  targets.  A  major  reason 
for  its  success  has  been  the  fact  that  the  EKF  has  offered  a  reasonable  compromise  between 
real-time  operation  and  accurate  performance  in  many  nonlinear  problems.  On  the  other  hand, 
the  EKF  is  a  suboptimal  and  completely  heuristic  algorithm  whose  efficiency  varies  from  case  to 
case.  For  instance,  the  EKF  is  unstable  in  situations  that  involve  acute  maneuvering,  missing 
measurement,  low  SNR,  multipath,  and  many  other  situations  where  the  posterior  distribution 
may  not  be  approximated  well  enough  by  a  Gaussian  distribution.  It  is  very  difficult,  if  even 
possible,  to  develop  rigorous  evaluation  metrics  for  assessing  the  quality  of  data  processing  based 
on  the  EKF  technology.  Improvements  to  the  EKF  (e.g.  iterated  EKF)  work  satisfactory  in  a 
number  of  important  applications  where  the  EKF  fails,  but  still,  it  is  difficult  to  overcome  the 
fundamental  limitations  of  the  EKF  algorithm  that  stem  from  its  dependence  on  the  assumption 
that  the  posterior  distribution  may  be  well  approximated  by  a  Gaussian  distribution.  The  typical 
posterior  density  built  upon  the  realistic  IR  image  is  shown  in  Figure  17,  Section  7  .  One  sees 
that  it  has  multi-peak  form  and  hence  is  very  far  from  being  Gaussian.  Our  experiments  show 
that  EKF  is  typically  fails  for  this  kind  of  data. 


2.1.4.  An  Alternative  TbD  Method  -  Optimal  Nonlinear  Filtering.  In  spite  of  the 
aforementioned  shortcomings,  the  above  outlined  methods  remain  the  basis  for  the  great  majority 
of  existing  signal  processing  systems.  In  particular,  until  recently  no  other  information  technology 
was  able  to  effectively  compete  with  the  EKF  in  target  tracking. 

However  the  situation  has  changed  with  recent  advances  in  the  mathematical  theory  and  algorith¬ 
mic  support  for  optimal  nonlinear  filtering  (ONF).  These  advances  coupled  with  improvements  to 
modern  digital  hardware  technology  make  optimal  nonlinear  filtering  an  attractive  alternative  to 
the  multi-dimensional  matched  filter,  dynamic  programming  based  algorithms  and  EKF  in  many 
practical  and  important  applications.  These  include  signal  processing  for  infrared  and  acoustic 
sensors,  imaging  radar  and  sonar  (including  SAR  and  SAS),  and  other  passive  and  active  sensors 
[17,  45,  49], 

Advanced  optimal  nonlinear  filters  can  now  provide: 

•  real-time  operation; 

•  optimal,  theoretically  sound  solutions  of  the  full  nonlinear  problem; 

•  distributional  versatility  (no  constraints  on  the  form  of  prior  or  posterior  probability  distri¬ 
butions); 
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•  superior  accuracy  and  robustness; 

•  facility  to  effectively  incorporate  realistic  physical  models; 

•  explicit  quantitative  performance  metrics  (exact  error  estimates,  confidence  areas,  etc.). 

Our  analysis  of  advanced  algorithms  based  on  ONF  technologies  shows  great  promise  in:  track- 
before-detect  of  unresolved  targets  in  low  S(N+C)R  (up  to  -6dB  after  preprocessing);  fusion 
of  imaging  and  kinematic  data  for  target  identification;  tracking  of  agile  extended  targets  in 
cluttered  environment  as  well  as  with  certain  other  applications. 

We  argue  that  just  as  the  EKF  superseded  the  a  -  /3  -  7  trackers,  so  the  optimal  nonlinear  filter 
is  set  to  replace  the  EKF  as  the  dominant  tracking  technology  within  10-15  years. 

In  Section  5  we  describe  the  ONF  techniques  that  will  be  used  in  defining  and  developing  the 
appropriate  technology  for  application  in  an  end-to-end  IRST  signal  processing  architecture. 
The  architecture  will  be  sufficiently  flexible  to  allow  for  an  adaptive  optimization  of  the  sig¬ 
nal/discrimination  processors  utilizing  existing  engineering  parameters.  This  adaptive  optimiza¬ 
tion  will  use  the  output  of  an  EO/IR  sensor  diagnostic  tool  to  utilize  current  meteorological  and 
environmental  information  as  well  as  current  intelligence  on  likely  threat  scenarios. 


3.  Signal  and  Observation  Modeling 

It  is  assumed  that  a  sensor  has  m-component  resolution  capability  (for  radar  typically  m  =  6,  for 
IR/EO  m  =  4).  By  r  =  (rq, . . . ,  rm)  will  be  denoted  a  phase  coordinate  vector  (for  IRST  typically 
angles  and  angle  velocities  of  an  object  in  a  certain  coordinate  system).  The  sensor  carries  out 
a  periodic  surveillance  of  definite  area  in  m-dimensional  domain  If  C  P,  where  Rm  is  the 
m-dimensional  Euclidean  space.  After  standard  preprocessing  the  result  of  one  observation  step 
is  a  frame  of  measurements, 

Z(n)  =  S(ri) +  b(ri) +ri(n)  n  =  1,2,...,  or 
Z(n)  =  \\Zi(n)\\  =  \\Si(n)+bi(n)  +  r)i(n)\\,  i  =  l,...,N, 

where  S(n)  =  ||Si(n)||  is  a  signal  from  target,  b(n)  =  ||6j(n)||  is  an  exterior  background  (clutter), 
and  rf  =  ||?/i(n)||  is  a  noise  of  the  sensor.  (Here  we  assume  that  after  preprocessing  sampling  of 
data  is  done  in  discrete  points  di,  i  =  1, . . .  ,  N,  uniformly  in  the  area  P71  (i  is  a  pixel).) 

The  noise  is  assumed  to  be  zero  mean  and  uncorrelated  in  both  time  n  and  space  i,  Er)i(n)  =  0, 
Era  («)  =  (E  is  a  symbol  of  expectation).  The  clutter  is  defined  as 

(3.2)  bi{n)  =  b(ri  +  S(n),n) 

where  b(r,  n )  is  a  function  describing  the  background  (spatial  distribution  of  the  clutter)  after 
preprocessing  in  the  point  r,  and  S(n)  =  (Sl(n), . . . ,  5m(n))  is  an  unknown  current  bias  of  sensor 
coordinate  system  with  respect  to  the  reference  one  (due  to  the  jitter). 
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The  signal  component  is  modeled  as 

k(n) 

(3.3)  Si(n)  =  Aj(n)h{ri  +  S(n)  -  r^n)) 

i=i 

where  h(r)  is  a  normalized  sensor  function;  k(n)  is  an  unknown  total  number  of  targets  at  the 
moment  n;  Aj(n)  and  r,(n )  are  unknown  signal  intensity  and  coordinates  of  the  jth  target, 
respectively. 

It  is  assumed  that  the  clutter  b(r,  n)  has  a  relatively  big  spatial  variance  (the  change  of  b(r,  n) 
between  two  nearest  values  of  r*  is  comparable  with  the  maximum  value  of  b(r,n)).  Besides, 
b(r,n)  varies  locally  as  fast  as  the  signal  function  does  (in  spatial  coordinates).  Even  in  cases 
where  these  assumptions  do  not  exactly  hold,  the  algorithms  that  rely  on  them  will  be  robust, 
which  is  the  most  important  requirement.  The  second  assumption  shows  that  if  only  spatial 
information  is  used,  then  the  background  may  be  interpreted  as  a  target.  Thus  spatial  filtering 
alone  is  not  sufficient  for  clutter  suppression  and  temporal  filtering  is  needed.  We  always  assume 
that  b(-)  is  an  arbitrary  unknown  function  of  r  and  slowly  changing  function  in  n:  there  exists 
such  T  that  \b(r,n  +  T)  -  b(r,n)\  <  av.  The  latter  assumption,  which  often  holds,  shows  that 
variations  of  the  clutter  in  time  are  caused  mainly  by  uncontrolled  vibrations  8{n)  of  the  sensor. 
These  vibrations  are  unknown  and  unpredictable  except  for  a  natural  restriction,  |6(n)|  <  A, 
imposed  on  the  absolute  value  of  bias  of  coordinate  system.  No  assumptions  on  statistical 
behavior  of  clutter  is  made.  It  is  our  belief  that  such  popular  models  as  homogeneous  random 
field,  especially  Gaussian,  are  valuable  only  for  purely  academic  research  and  lead  to  highly 
non-robust  filtering  algorithms. 

Perhaps  the  most  important  feature  of  our  approach  to  algorithm  design  is  that  we  refuse  to  use 
any  artificial  and  unreliable  statistical  models  of  b(r,n),  8(n),  k(n),  and  Aj(n).  The  algorithms 
developed  on  the  basis  of  such  models  fail  even  for  small  deviation  of  a  model  from  reality. 
The  essence  of  the  new  suggested  approach  is  the  development  of  algorithms  that  are  invariant 
and/or  adaptive  with  respect  to  prior  uncertainty.  The  specific  feature  of  the  problem  is  its 
extremely  high  dimensionality.  The  value  of  N  can  be  of  the  order  of  106  -  108  and  the  total 
number  of  targets  can  be  several  hundreds.  Thus,  along  with  the  mathematical  issue  of  data 
processing,  serious  attention  should  be  paid  to  the  computational  complexity,  parallelization  and 
HPC  realization. 


4.  Clutter  Suppression 

Two  classes  of  algorithms  for  background  filtering  are  proposed.  In  either  case  there  are  two 
basic  problems  to  be  solved:  (a)  to  transform  the  sequence  of  input  frames  into  the  new  frame 
such  that  the  clutter  would  be  reduced  and  the  signal  would  be  preserved  (clutter  removal); 
(b)  for  every  n  the  location  of  the  coordinate  system  of  the  sensor  (i.e.  the  bias  8(n))  should 
be  estimated  with  maximum  possible  accuracy  (jitter  compensation).  Below  we  propose  two 
algorithms  that  allow  the  first  problem  to  be  solved  effectively.  Jitter  compensation  algorithms 
will  be  developed  in  the  near  future. 
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4.1.  Nonparametric  Method.  The  first  class  is  relied  on  the  nonparametric  regression  ap¬ 
proach  to  estimation  of  the  function  b(r,  n ).  Our  goal  is  to  build  a  nonparametric  clutter  estimate 
such  that  the  residuals  between  the  original  data  and  its  smoothed  version  (estimate)  would  be 
reasonably  approximated  by  signal  plus  noise  models.  That  is,  the  estimate  b(r,  n )  of  the  b(r,  n) 
should  be  built  in  such  a  way  that  in  the  filtered  frame  Z(r,n )  =  Z(r,n )  -  b(r,n )  the  signal 
would  be  preserved,  while  the  clutter  would  be  removed  almost  completely. 

Kernel  methods  provide  a  powerful  tool  for  such  analysis  due  to  both  theoretic  optimality  (see 
[19,  23])  and  computational  transparency.  In  addition  these  methods  are  invariant  to  statis¬ 
tical  properties  and  variations  of  clutter.  Kernel  estimators  are  weighted  moving  averages  of 
observations 

S(r'n)  =  rrnv  £  *(»•.»)*  (V--^)' 

where  Ni  are  window  sizes  in  corresponding  directions,  K(-)  is  a  deterministic  m— dimensional 
kernel  (recall  that  r  =  (ri, . . . ,  rm)).  A  product  kernel  with  univariate  kernels  Ki,  i  =  1, . . . ,  m, 
provides  a  popular  example. 

Note  that  in  the  2D  case  (m  =  2)  a  so  called  ‘local  mean  removal’  procedure  used  for  clutter 
rejection  in  a  series  of  papers  by  Reed  et  al.  [38,  39]  is  equivalent  to  a  kernel  estimator  with 
a  product  kernel  K  generated  by  uniform  kernels  K\{x)  =  K<2(x)  =  0.5/{|x|<i}.  It  is  shown 
in  Section  4.1  that  the  proposed  adaptive  nonparametric  methods  based  on  recent  advances  in 
nonparametric  estimation  (see  e.g.  [18],  [27])  perform  much  better  than  the  aforementioned  ‘local 
mean  removal’  method. 


4.2.  Semiparametric  Filtering.  The  second  method  (semiparametric)  is  based  on  the  adap¬ 
tive  spatio-temporal  auto-regression 

n  L 

(4.1)  Zi(n)  =  Zi(n)-  ^aj(r)Z«(n) 

r—n-T  1—0 

where  Z(n)  =  ||£,(n)||  is  the  filtered  frame  of  input  data  Z{n)  =  ||Zj(n)||;  Zu(n)  are  the  values 
of  Zi(n)  in  some  vicinity  (spatial  filtering  window)  of  points  r*;  L  is  the  number  of  points  rj  in 
the  window;  ai(r)  are  some  coefficients. 

The  major  problem  is  the  choice  of  optimal  coefficients  o;(r).  They  must  be  calculated  adap¬ 
tively  to  guarantee  minimum  of  empirical  mean-square  (MS)  value  of  the  filtering  residual  noise 
(internal  noise  and  background  residual)  for  every  time  moment  n.  Such  a  criterion  provides 
minimum  of  MS  residual  noise  in  the  frame  for  every  time  n  and  is  invariant  with  respect  to 
statistical  properties  of  the  background  b(ri  +  <5(n),n)  and  noise.  The  algorithm  consists  of: 

1.  Computation  of  empirical  correlation  matrix  (TL  xTL)  for  spatial  window  L  and  temporal 
window  T. 

2.  Calculation  of  the  optimal  weights  aj(r). 
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3.  Linear  transform  of  Zi{n)  into  Zi(n)  according  to  (4.1). 

In  general  the  algorithm  is  nonlinear  because  of  the  nonlinear  dependence  of  at(r)  on  Zi.  This 
algorithm  allows  us  to  suppress  any  background  regardless  of  its  variation  in  time.  However, 
it  solves  only  the  first  part  of  the  problem  -  it  does  not  estimate  a  drift  of  sensor  coordinate 
system  to  correct  coordinate  measurements  on  the  stage  of  detection  (jitter  compensation).  An 
algorithm  that  allows  us  to  estimate  6(n)  and  hence  to  compensate  jitter  will  be  developed  in 
the  near  future. 


5.  Optimal  Nonlinear  Filtering  for  TbD 

5.1.  Modeling  for  TbD.  Our  philosophy  in  addressing  the  low  S(N+C)R  detection  problem 
(regardless  of  the  type  of  a  sensor  being  used)  is  that  the  farther  the  thresholding  operation 
may  be  postponed  the  better.  In  the  context  of  track  before  detect,  we  assume  that  the  data 
are  collected  at  discrete  times  while  the  target  dynamics  is  a  continuous  time  process.  To  be 
specific,  assume  that  the  measurements  Z(k,  r )  =  Z(tk,  r )  are  collected  at  discrete  time  moments 
tk,  k  =  0,1,2,...  and  the  relationship  between  the  observation  and  target  location  is  modeled 
by  a  nonlinear  measurement  equation  of  the  form 

(5.1)  Z(k,  r )  =  S(k,  Xk,  r)  +  b{k,  r )  +  V(k,  r ) , 

where,  as  before,  r  represents  the  spatial  coordinate  in  the  phase  space,  b(k,  r)  is  a  sequence  of 
deterministic  (unknown)  functions  representing  cluttered  background,  by  Xk  is  denoted  the  true 
location  of  the  target  at  time  t  —  tk  and  by  V(k,r )  is  denoted  the  measurement  noise  process 
(sensor  noise).  For  simplicity  we  first  suppose  that  there  may  be  only  one  target  in  the  scene,  i.e. 
the  signal  S(k,Xk,r)  is  described  by  (3.3)  where  k(n)  =  1.  Also  we  will  neglect  the  platform 
instability  (jitter)  assuming  S(n)  =  01.  Under  these  conditions 

(5.2)  S(k,Xk,r)  =  A(k)h(r-Xk), 

where  h(r )  is  a  normalized  sensor  function  and  A(k)  is  an  unknown  signal  intensity. 

We  further  define  Zk(r )  =  (Z(l,  r), . . . ,  Z(k,  r))  to  be  the  concatenation  of  all  measurements  up 
to  time  4  in  the  space  point  r. 

The  expected  range  of  possible  “behaviors”  of  the  target  is  modeled  as  a  Markov  process.  Often 
this  process  may  be  well  described  by  a  randomly  perturbed  multi-dimensional  linear  or  nonlinear 
dynamical  system 

Xt  =  f(Xt)  +  aWt,  t>  0, 

where  Wt  is  another  noise  process  that  describes  uncertain  and  unpredictable  target  motion;  /(•) 
is  a  known  function.  Modeling  of  the  state  and  observation  is  normally  based  upon  physics.  The 
models  for  the  state  process  usually  represent  the  a  priori  knowledge  about  physical,  tactical, 
etc.  characteristics  of  the  target  while  the  observation  model  is  based  upon  the  physics  of  sensors, 
clutter  structure,  operational  environment,  and  so  forth. 


1  Tli is  is  the  case  if  the  instability  is  compensated  either  by  electro-mechanical  stabilizators  or  by  estimating. 
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5.2.  Optimal  Spatio-Temporal  Nonlinear  Filtering:  The  Basic  Algorithm.  The  output 
of  the  ONF  is  the  functional  time  series  {7rfc(r)}fc>i  of  joint  posterior  densities  (JPD)  which 
measure  the  likelihood  at  the  time  moments  4  that  the  vector  of  target  features  parametrized 
by  Xk  is  in  close  proximity  of  the  grid  point  r. 

Computation  of  the  JPD  usually  splits  into  two  separate  procedures:  on  every  time  step,  spatio- 
temporal  filtering  for  clutter  removal  is  done  first  (see  Section  4),  and  then  the  spatio-temporal 
nonlinear  filtering  is  performed  to  estimate  the  target  location.  After  clutter  suppression  proce¬ 
dure  the  clutter  is  decorrelated.  We  thus  can  include  it  into  the  noise  component  V(k,r)  (see 
(5.1))  and  concentrate  only  on  the  second  procedure. 


It  is  a  standard  fact  (see  e.g.  [24],  [42])  that  the  JPD  is  given  by  the  formula 

Pk(r) 


(5.3) 


7 Tfc(r)  = 


Jpk(r)dr  ’ 

where  the  function  pk(r),  called  the  unnormalized  filtering  density  (UFD),  propagates  in  time 
according  to  the  following  recursive  equation  of  the  predictor-corrector  type: 

(5.4)  pk(r)  =  exp{(fi-1/i(r,  •),  R~lZ{k,  •))  -  \  \ R~lh{r,  -)| 2}Tkpk.x{r),  k  =  1,2,..., 

where  R  is  the  covariance  of  the  observation  noise  V,  and  the  predictor  Tkpk_x(r)  is  a  solution  of 
the  Fokker-Planck-Kolmogorov  equation  corresponding  to  the  state  process  subject  to  the  initial 
condition  pfc_i(r),  i.e.  u(t,r)  =  Ttpk-i(r)  is  a  solution  of  the  equation 


(5.5) 


du(t,  r ) 
dt  - 

*j=i 

u(0,r)  =pfc_i(r), 


=  I  “  £  ^(/i(r* Ht,r)),  t  e  (4,4- 1], 


d 


dr$r3 


i=l 


dr i 


where  {ay}  is  the  covariance  of  the  state  noise  and  m  is  the  dimensionality  of  Xt.  The  posterior 
moments  of  Xk  =  Xtk  can  be  computed  now  by  integration  of  the  respective  polynomials  against 
the  JPD  7rfc(r).  For  example,  the  best  mean  square  estimate  of  Xk, 


dr. 


Recall  that  we  consider  a  general  problem  assuming  that  r  belongs  to  a  rn- dimensional  Euclidean 
space  (phase  space)  (our  main  concern  in  IRST  problem  is  m  =  4  if  the  problem  is  solved  in 
the  initial  “resolution”  space  (angles  and  angular  velocities)  or  m  =  2  if  trajectories  are  projected 
on  the  plain). 

Now,  let  {et{r)}i£N  be  a  complete  orthonormal  system  in  L2 (Rm ) .  Projecting  (5.1)  on  this  basis, 
we  can  rewrite  the  observation  process  in  the  coordinate  form: 

(5.6)  Zi  =  hl(Xtk)  +  Vk\  *  =  1,2,...,  *  =  1,2,...  ,K, 

where 

Zek=  [  Zk(r)ee{r)dr,  he{y)  =  [  h{y,r)ee{r)  dr,  Vk  =  f  Vk{r)et{r)dr 

JUrn  Jl&m 

and  Vk  are  independent  Gaussian  random  variables. 
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Let  F(-)  be  a  given  function  such  that  E  |F(Xt)|2  <  oo,  Vi  >  0.  In  general,  the  optimal  nonlinear 
filtering  deals  with  obtaining  the  best  mean  square  estimate  of  F(Xtk),  given  the  measurements 
Zk[r )  —  (. Zi(r),...,Zk(r )),  r  e  Rm.  This  estimate  is  called  the  optimal  filter  and  will  be 
denoted  F(k). 


By  Tt  denote  the  solution  operator  for  the  Fokker-Planck-Kolmogorov  equation  corresponding 
to  the  process  X.  The  optimal  filter  F(k)  =  E[F(Xtk)  \  Zk]  is  given  by  the  formula 


(5.7) 


_  Lr>P*(r)F(r)dr 
()~  L.Pt(r)dr  • 


where  the  UFD  pk(r)  was  introduced  above.  This  density  obeys  the  recursive  relation 


(5.8) 


p0(r)  =  tt  (r), 

f  OO  A  00  . 

Pk(r)  =  exp  A  Y,  h‘(r)Z‘  -  §  Y.  I^MI2  Tzto-d  r),  k  =  l,2.... 


1=  1 


t=l 


Additional  notation: 

Qmn  ~  (emiTACn)j  Qmn  n)>  Qmn  ^  ®mjTAfin)> 

'0m(  0)  =  (P>®m))  Fm  —  (F,em), 

where  {F,  g)  =  fRm  F(r)g(r)  dr.  Note  that  these  quantities  can  be  computed  off-line  before  any 
observations  become  available. 

Let  L  =  {£:  £  <  M}  and  Fk  be  the  approximation  to  the  optimal  filter  F(k). 

We  propose  the  following  algorithm  for  computing  the  optimal  nonlinear  filter  F(k),  which  will 
be  called  the  Spectral  Separation  Scheme  or  the  S 3  algorithm. 

1.  Set  a  cut-off  level  M  for  the  number  of  basis  elements. 

2.  Before  the  observations  become  available  compute: 

Po(r)  =  M®)ee(r)  and  %[f]  =  ^  ^(0 )Ft. 

t<M  t<M 

3.  When  the  measurements  ( ),  k  =  1,2,...  ,K,  £  e  L  C  L  become  available: 

a) :  compute  recursively 

V’m(O)  =  V’m(O),  ^(k)=Y^Qmn(Zk)fPn(k-l),  k=  1,2,...; 

n<M 

0™„(Z^)=?m»  +  Aj]ZiE9L  +  A2  E  +  A-iyi; 

t£L  ift'AjZL  leL 

b) :  compute 

ph(r)  =  E  **[F]  =  E  Fk  =  **[F]/f  *[1]. 

1<M  1<M 
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The  block  diagram  of  the  Spectral  Separation  Scheme  (S'3  algorithm)  is  shown  in  Figure  2.  It 
is  clear  that  the  above  algorithm  allows  for  an  adaptive,  dynamic  interplay  between  data  and 
computations.  Indeed,  it  does  not  involve  solving  of  PDE’s  on-line  and  is  recursive  in  time. 
Furthermore,  it  is  also  spatially  recursive.  These  properties  are  very  important  since  they  allow 
adaptive  sequential  multi-resolution  filtering  to  be  performed.  Observe  that  in  order  to  achieve 
the  effect  of  multi-resolution  one  has  to  use  local  hierarchical  bases  {e;},  for  example  wavelet 
bases  [25]. 


Yik-l) 

yr 

observation 

Optimal  Nonlinear  Filter 

V,(fc)=X  QnV (*-l) 

n<M 

Q 

W) 

Estimation 

Pic(r)  =  '2dV/i(k)el(r) 

1<.M 

W 

Q  -  precomputed  matrix 


Figure  2 .  S 3  algorithm  (on-line  part) 


5.3.  Multi-Dimensional  Spatio-Temporal  Matched  Filter  as  a  Special  Case  of  Nonlin¬ 
ear  Filter.  For  the  sake  of  simplicity  let  us  consider  the  basic  formula  for  updating  the  filtering 
density  in  the  optimal  nonlinear  filter  in  the  3— D  case  where  the  target  motion  is  projected  on 
the  plain.  The  modification  for  the  (m  +  1)-D  case  with  arbitrary  m,  particularly  for  m  =  4, 
is  straightforward.  The  grid  space  will  be  denoted  by  {xij}.  Suppose  that  the  solution  of  the 
Fokker-Plank-Kolmogorov  equation  (the  mean  motion  dynamics  of  the  target)  is  given  by  the 
operator  Ta  and  observations  are  made  at  time  moments  4  =  on  the  observation  space  grid 
{xij}  with  additive  Gaussian  white  noise: 

(5.9)  Z^h)  =  h(Xtk  —  x^ )  +  <x£ij{tk)i 

where  h(-)  is  the  target  signature  signal  and  £^(4)  are  i.i.d.  standard  Gaussian  variables. 

Then  by  (5.4)  the  updating  formula  is 

Pk{x)  —  ^  %ij)Zij{tk)  —  2^2  ^  >  h(x  —  Xij)  |TAPfc-i(:c)  . 

i,j  hi 
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In  the  above  formula,  the  part  TApk-i(x)  is  the  prediction  part,  i.e.  the  best  possible  guess  (based 
on  known  target  trajectory  dynamics),  what  the  probabilities  to  find  the  target  in  various  areas 
of  the  observation  region  would  be  if  there  were  no  observations  available  at  the  time  moment 
tk.  The  exponent  part  is  the  correction  part,  i.e.  after  receiving  the  observation  at  time  tk,  the 
prediction  part  is  amplified  or  attenuated,  depending  on  the  expression  inside  the  exponent. 

One  notable  fact  is  that  the  amplification  or  attenuation  is  directly  related  to  the  output  of  the 
well-known  in  engineering  community  single  frame  matched  filter  (spatial  (2— D)  matched  filter) 

mk{x)  =  h(x  -  Xij)Zij(tk). 
i,j 

Indeed,  the  filtering  densities  are  amplified  only  at  those  points  where  the  output  of  the  matched 
filter  exceeds  the  value  of 

|E2  =  \ 

i,3 

and  attenuated  elsewhere.  Here  E2  is  the  energy  of  the  signal  which  of  course  does  not  depend 
on  x. 


This  establishes  a  strong  connection  between  a  spatial  matched  filter  and  the  optimal  nonlinear 
filter.  Even  more  importantly,  the  optimal  nonlinear  filter  is  equivalent  to  the  spatio-temporal 
matched  filter  (or  assumed  velocity  matched  filter)  if  the  movement  of  the  target  occurs  along  a 
deterministic  trajectory.  In  that  case  the  operator  TA  degenerates  into  an  operator  that  shifts 
the  values  of  the  function  along  the  known  deterministic  trajectory.  For  simplicity,  let  us  assume 
that  the  target  moves  in  a  straight  line  with  known  constant  velocity  v.  Then 

TAp(x)  =p{x-v  A), 


and  it  is  easy  to  see  that  for  any  k  >  1 
fc-i 


1  fc_1  1  *-1 

pk{x)  =  ck  exp  |  xij  -  ivA  )Zij(tk)  -  ^  Y2  h(x  ~  Xii  ~  foA)/Po(z  ~  kvA). 


e=o  i,j 


e-o  i,j 


Notice  that  the  first  part  in  the  exponent  in  the  above  formula, 

fc-i 

Mk(x)  =  EE  h ( x  x^  £vA)Zij(tk) , 

e=o  i,j 

exactly  coincides  with  the  usual  assumed  velocity  matched  filter  and  the  second  part 

2^EE^- “  ^A)2  =  ^2k 

t=0  i,j 

is  nothing  but  the  accumulated  SNR.  Thus,  if  we  assume  that  the  initial  density  is  a  delta- 
function,  then  statistical  inference  base  on  Mk(x)  and  pk(x)  will  be  the  same.  Particularly,  the 
points  of  maxima  of  Mk(x )  and  pk(x)  coincide. 

A  similar  argument  can  be  applied  to  any  dimensions.  (Recall  that  in  general  we  are  interested 
in  the  (m  +  1)-D  case,  r  =  (rx, . . .  ,rm),  (m  +  l)-th  component  is  time.)  The  final  result  is 
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exactly  the  same:  in  case  of  deterministic  trajectory  the  (m+l)—D  matched  filter  and  the  optimal 
spatio-temporal  nonlinear  filter  coincide  with  the  accuracy  to  the  deterministic  transformation. 

We  stress  that  the  developed  optimal  spatio-temporal  nonlinear  filter  does  more 
than  the  spatio-temporal  matched  filter  does.  It  is  a  generalization  of  the  multi¬ 
dimensional  spatio-temporal  matched  filter  and  coincides  with  it  for  targets  that 
move  along  deterministic  trajectories.  In  general ,  however,  the  ONF  predicts  the 
trajectory  and  allows  us  to  align  frames  coherently  even  for  acutely  maneuvering 
targets  when  the  (m  +  1  )-D  matched  filter  fails. 

6.  Track  Appearance/Disappearance  Detection 

6.1.  Preliminaries.  A  problem  of  detecting  target  tracks  that  occur  (appear  and  disappear)  at 
a  priori  unknown  points  in  time  is  a  typical  abrupt  change  detection  problem.  The  change-point 
detection  problem  has  been  actively  researched  during  the  last  three  decades  (see,  e.g.,  [22,  34, 
37,  41,  46,  44,  52,  53,  54]  and  references  therein).  Heuristic  procedures  such  as  Shewhart’s  charts 
and  some  modifications  appeared  in  the  late  twenties  and  early  thirties.  Detection  procedures, 
which  are  in  current  use,  were  initiated  by  Girshik  and  Rubin  [22]  and  Page  [34]  for  the  problem  of 
detecting  a  change  in  a  mean  of  i.i.d.  Gaussian  sequence  and  by  Shiryaev  [44]  in  the  general  (but 

1.1. d.)  case.  There  are  two  major  competitive  procedures:  the  Shiryaev-Roberts-Girshik-Rubin 
algorithm  [41,  37,  54,  44]  and  the  Page’s  (or  CUSUM  -  cumulative  sum)  procedure  [34,  46,  55]2. 

There  are  two  major  disadvantages  to  both  algorithms: 

1.  Only  one  change  in  data  is  assumed,  i.e.  in  our  context  the  moment  of  target  disappearance 
is  ignored. 

2.  The  aforementioned  change-point  detection  procedures  have  optimal  properties  only  in  the 

i.i.d.  case  in  the  following  specific  sense:  they  minimize  the  mean  detection  delay  under 
constraints  on  the  mean  time  between  false  alarms  for  exactly  specified  pre-change  and 
post-change  distributions. 

These  drawbacks  make  the  conventional  change-point  detection  procedures  impractical  in  multi¬ 
target  surveillance  problems  that  involve  non-i.i.d.  observations  and  when  one  needs  not  only  to 
detect  a  particular  target  (with  unknown  position)  but  also  to  discriminate  between  false  alarms 
and  true  tracks.  One  appropriate  criterion  for  our  purposes  is  to  fix  probabilities  of  false  alarm 
and  detection  in  a  fixed  size  sliding  window  and  to  minimize  the  detection  delay.  Or,  alterna¬ 
tively,  one  may  require  to  detect  the  track  with  fixed  or  maximum  probability  during  the  fixed 
time  and  with  minimum  detection  delay.  Also  it  is  imperative  to  design  the  decision  statistics 
that  will  account  for  unknown  target  location  and  moment  of  its  disappearance.  Otherwise  the 
problem  is  meaningless.  Neither  of  the  above  mentioned  algorithms  may  be  directly  applied  in 
this  context.  However,  we  will  show  that  the  CUSUM-type  procedure  and  the  quasi-Bayesian 
(Shiryaev-Roberts-Girshik-Rubin  type)  procedure  with  specially  designed  thresholds  are  useful 


2See  also  [40,  54]  for  reasonable  modifications  and  ideas  that  are  useful  for  surveillance  systems. 
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tools,  especially  in  the  problem  of  detection  of  a  single  target  that  appears  and  disappears  at  un¬ 
known  times.  Moreover,  it  turns  out  that  the  thresholds  can  be  chosen  such  that  the  predefined 
false  alarm  rate  is  fixed  for  a  large  class  of  noise  and  signal  models,  not  necessarily  for  the  i.i.d. 
models.  This  innovation  is  very  important  for  our  applications,  since  it  removes  the  restrictive 
i.i.d.  assumption,  which  is  typical  for  previous  work  in  change-point  detection.  In  particular,  it 
allows  us  to  design  thresholds  for  adaptive  algorithms. 

Now  we  summarize  the  basic  requirements  that  should  be  satisfied  in  realistic  problems  of  de¬ 
tection  of  target’s  track  appearance/disappearance: 

•  detection  algorithms  should  be  invariant  relative  to  the  completely  unknown  moments  of 
abrupt  change  (the  prior  distributions  of  appearance/disappearance  moments  are  unknown); 

•  detection  algorithms  should  be  adaptive  and  use  estimates  of  unknown  target  location  based 
on  TbD; 

•  the  frequency  of  false  detections  or  the  false  alarm  probability  is  limited  at  the  specified 
level; 

•  the  probability  of  correct  detection  should  be  maximal  during  the  fixed  time  interval; 

•  the  mean  detection  delay  should  be  made  as  small  as  possible. 

We  propose  the  sequential  algorithms  that  meet  all  these  requirements  (and  even  more,  see  Algo¬ 
rithm  3  in  Section  6.5).  The  algorithms  are  based  either  on  the  generalized  CUSUM  statistic  or 
on  the  quasi-Bayesian  statistic.  These  statistics  are,  in  essence,  score  functions  that  characterize 
the  likelihood  ratio  of  hypotheses  on  signal  presence  and  absence.  Specifically,  the  generalized 
CUSUM  statistic  is  the  likelihood  ratio  maximized  over  unknown  moments  of  appearance  and  dis¬ 
appearance  and  the  latter  statistic  is  the  average  likelihood  ratio  with  respect  to  flat  (improper) 
prior  distributions  of  these  moments.  The  statistics  are  computed  either  in  a  sliding  window  of  a 
fixed  size  or  on  the  semiinfinite  interval  (up  to  the  current  moment).  In  both  cases  the  decision 
statistics  are  compared  with  thresholds  that  depend  on  a  given  false  alarm  rate  (probability  of 
false  alarm  or  frequency  of  false  alarms).  It  turns  out  that  in  many  cases  the  generalized  CUSUM 
procedure  may  be  reduced  to  the  conventional  CUSUM  without  loss  of  statistical  performance. 
To  be  precise,  if  the  major  goal  is  to  detect  the  fact  of  target’s  appearance  (but  not  the  fact  of 
its  disappearance),  than  the  unknown  moment  of  target’s  disappearance  should  be  considered 
as  a  nuisance  parameter.  In  this  case  the  CUSUM  is  the  optimal  algorithm,  i.e.  the  unknown 
moment  of  target’s  disappearance  does  not  affect  the  structure  of  the  algorithm,  but  only  its 
performance  (see  Section  6.3  and  Section  6.4  for  details).  The  same  is  true  for  the  quasi-Bayes 
procedure.  If,  however,  one  has  to  detect  multiple  targets,  then  it  is  important  to  detect  both 
track  appearance  and  track  disappearance  (as  soon  as  possible).  Moreover,  in  this  case  it  is  also 
important  to  have  a  special  logic  to  discriminate  between  false  alarms  and  true  tracks.  In  this 
“multi-target  context”  the  conventional  CUSUM- type  procedures  do  not  solve  the  problem  but 
still  may  be  used  as  a  part  of  more  sophisticated  algorithms  (see  Section  6.5  for  some  discussion). 

It  should  be  mentioned  that  the  problem  of  signal  detection  with  random  appearance  and  dis¬ 
appearance  times  in  ideal  conditions  (complete  prior  information)  was  solved  by  Tartakovsky 
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[56].  The  results  of  this  work  show  that  the  optimal  Bayes  solution  requires,  in  general,  com¬ 
putation  of  the  5-dimensional  sufficient  statistic  and  is  completely  impractical.  Thus  even  the 
availability  of  the  required  prior  information,  particularly  knowledge  of  prior  distributions  of 
appearance/disappearance  moments  (which  are  usually  unknown),  does  not  help  in  practice. 


6.2.  Problem  Formulation  and  Model  Description.  Let  Yk,  k  —  1, . . .  ,n,  denote  the  data 
obtained  up  to  the  current  time  n.  Generally,  the  data  Yk  =  {Zk(ri),  i  =  1, . . . ,  N}  represent  the 
whole  “spatial  set”  collected  from  the  frame  of  observations  (after  preprocessing)  at  time  k  or 
the  part  of  these  data  (in  some  “channels”).  If  the  clutter  removal  procedure  is  first  performed, 
then  (after  spatio-temporal  filtering)  the  clutter  is  decorrelated.  Thus,  it  is  reasonable  to  assume 
that  {Yk}  are  i.i.d.  with  the  density  p\{y)  if  the  target  is  present  and  with  the  density  p0(y)  if  it 
is  absent3.  The  density  p\{Yk  |  9k)  depends  on  (in  general  unknown)  vector  parameter  9k,  which 
characterizes  the  spatial  location  of  the  target  at  moment  k.  First,  for  the  sake  of  simplicity 
we  assume  that  the  target  dynamics  is  known  and  hence  the  9k  is  fixed  and  known.  The  exact 
models  with  the  account  of  TbD  and  corresponding  adaptive  algorithms  will  be  defined  later  on 
in  Section  6.8  and  Section  6.9  (see  also  Section  6.6  in  the  general  case). 

By  Ak  =  log  denote  the  log-likelihood  ratio  (LLR)  for  the  single  observation  Yk  and  by  A 
and  7  the  unknown  moments  of  target  appearance  and  disappearance,  respectively.  Under  our 
assumptions  the  data  Yi, . . . ,  Y\-i  are  i.i.d.  according  to  the  density  po(y),  Yx, . . . ,  F7  are  i.i.d. 
according  to  pi(y)  and  F7+i,  F7+2, . . . ,  n  again  follow  p0(y).  We  assume  nothing  about  A  and  7 
except  for  the  natural  constraint  7  >  A.  In  the  language  of  hypotheses  testing  the  problem  of 
the  track  detection  may  be  formulated  as  testing  of  the  hypotheses 

:  track  appears  at  A  =  l  and  disappears  at  7  =  m,  1  <  l  <  n,  m  >  l  +  1”; 

“ Ho  :  track  does  not  appear,  i.e.  A  ^  [l,n]”. 

Alternatively,  the  hypotheses  may  be  written  in  the  form 

“iru,7 :p(i?) =p^(y?) = n»,(n) x n*.«) x  n  • 

k~l  k=\  fc=7+l 

k=l 


6.2.1.  Generalized  Likelihood  Ratio  Statistic.  By  Ln  denote  the  generalized  LLR, 


Ln  :=  max  log 

a, 7 


Pop?) 


7 


max 

l<A<n 


7>A+1  k=X 


3The  derived  algorithms  can  be  easily  modified  and  generalized  for  more  general  models  that  include  correlated 
and  non-homogeneous  observations,  see  comments  below  and  Section  6.7. 
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It  is  easily  shown  that  the  statistic  Lk  satisfies  the  following  system  of  recursive  relations: 

(6.1)  Lk  =  max(Lfc_i,  Uk)  =  Lfc_i  +  ( Uk  -  Lk- 1)+, 

(6.2)  Uk  =  Afc  -I-  Uk_i,  k  =  1, . . . , 7i, 


with  the  initial  conditions  L0  =  U0  =  0.  Here  y+  denotes  a  non-negative  part  of  y,  i.e.  y+  = 
max(0,y).  Indeed, 


Ln  =  max 

7>1 


(  max  VAfc  =  max(f/i, ...,Un)=  max(Ln_i,  t/„) , 

J 


where  for  any  k  >  1 


k 

Uk  =:  max  V  As  =  max 

l<A<fc  ^ 
s= X 


At,  max 

l<A<fc-l 


fc-i 

Afc  +  As 

s=X 


=  A*  +  max(0,  Uk-i) . 


6.2.2.  Generalized  Average  Likelihood  Ratio  Statistic.  Introduce  the  following  statistic 


Gn  :=  max 

7 


I?^rfA=^§nexp(At)' 


It  is  easy  to  see  that 


Gn  —  max(jF?i,  R 2, . . . ,  Rn)  —  max((jn_i,  Rn)i  Gq  0, 

where  Rn  =  £)"=1  IIfe=A  eA*  ■  The  statistic  Rn  may  be  interpreted  as  the  likelihood  ratio  averaged 
over  the  flat  (uniform)  distribution  of  the  moment  of  track  appearance. 


It  is  easy  to  show  that  {Rn}  obeys  the  recursive  relation 
(6.3)  Rn  —  exp(An)(l  +  i?n— 1),  Ro  =  0  . 


Both  statistics,  Ln  and  Gn,  can  be  used  for  testing  the  above  hypotheses  either  in  the  fixed 
interval  or  sequentially. 


6.3.  Detection  Algorithm  1:  Fixed  Sliding  Window.  The  hypotheses  #1^,7  and  H0  are 
tested  in  a  sliding  window  of  length  T  at  each  current  time  instant  n,  i.e. 

A— 1  7  n 

“■Hi, A, 7  '■  PO^n-T+l)  =  PX.'rO^n-T+l)  =  PoO^k)  X  JJpi(Ffc)  X  jQ  Po(Tfc)”  , 

k=n~T+l  k= X  k= 7+I 

n 

“■ffo:pre_T+i)=poK-T+i)=  n  Mnr- 

k=n—T+l 

In  this  case  the  statistics  LT,n  and  Gx,n  are  determined  by 

(6.4)  LT<k  =  max^fc-x,  Uk), 

(6.5)  Uk  =  Afc  +  U£_i,  k  =  n  —  T+l,...,n, 
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and 

(6.6)  GT,k  =  max(Gr,fe_i,  Rk ) , 

(6.7)  Rk  =  eAfc(l  +  -Rfc-i))  k  =  n  —  T  +  l,...,n, 

with  the  initial  conditions  Lx,n-T  =  Un-T  —  Gt,u-t  =  Rn-r  =  0. 


According  to  the  adaptive  Bayesian  approach,  the  procedures  of  testing  the  hypothesis  Hi  against 
Hq  have  the  form 


(6.8) 


if  LT,n  >  a, 

otherwise, 


(6.9) 


f  1  if  GTtn  >  B, 
I  0  otherwise, 


where  a,  B  are  thresholds  that  are  defined  on  the  basis  of  the  given  probability  of  false  alarm 
and  d  =  1  stands  for  the  decision  on  target  presence  while  d  =  0  is  the  decision  on  its  absence. 
It  should  be  emphasized  that  in  this  subsection  we  pursue  the  goal  of  detection  of  the  target 
track  appearance  regardless  of  the  possibility  of  its  disappearance  in  the  analyzed  observation 
interval.  In  other  words,  the  unknown  moment  of  target  disappearance  is  considered  as  a  nuisance 
parameter.  A  different  setting  when  both  target’s  appearance  and  disappearance  should  be 
detected  will  be  considered  in  the  next  two  sections. 


We  now  notice  the  following  remarkable  property  of  the  tests  (6.8)  and  (6.9).  Since 


LT,n=  max  Uk, 

n— T+l<k<n 


we  have 

{d^n  =  0}  =  {Uk  <  a  for  all  n  —  T  +  1  <  k  <  n}  . 

This  shows  that  the  non-sequential  test  (6.8)  is  equivalent  to  the  sequential  test  which  is  defined 
by  the  stopping  time 

(6.10)  r0(T,  n)  =  min{n  -  T  +  1  <  k  <  n  :  Uk  >  n},  ra(T,  n)  =  oo  if  no  such  k. 

If  ra(T,  n)  <  oo,  then  the  hypothesis  Hi  is  accepted  (target  track  is  present),  while  if  r0(T,n)  = 
oo,  the  hypothesis  H0  (target  is  absent)  is  accepted.  The  same  argument  is  applied  to  show  that 
the  non-sequential  test  (6.9)  is  equivalent  to  the  sequential  one 

(6.11)  fjg(T,  n)  =  min{n  -T+l<k<n:Rk>  B},  tb(T,  n)  =  oo  if  no  such  k. 


The  sequential  procedures  (6.10)  and  (6.11)  allow  us  to  achieve  exactly  the  same  probabilities 
of  false  alarm  and  target  missing  as  the  initial  tests  (6.8)  and  (6.9),  respectively.  In  addition, 
they  have  an  advantage:  the  time  delay  in  signal  detection  is  less.  See  Figure  3  for  graphical 
explanation  of  the  above  argument. 
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Figure  3.  Detection  in  the  sliding  window:  the  statistic  Un  exceeds  the  threshold 
at  random  moment  ra,  which  is  less  than  the  size  of  the  sliding  window  T. 

We  also  observe  the  following  interesting  fact.  The  statistic  Uk  is  nothing  but  the  CUSUM 
statistic,  which  is  defined  in  the  standard  CUSUM  procedure  as 

k 

Uk  =  max  V]  As, 

A< k  ' 

S=\ 

i.e.  as  the  maximum  of  the  LLR  over  unknown  moment  of  signal  appearance  only.  In  turn,  the 
statistic  Rn  is  nothing  but  the  quasi-Bayesian  statistic,  which  is  the  limit  form  of  the  Bayesian 
statistic  (the  Shiryaev-Roberts-Girshik-Rubin  statistic).  This  shows  that  the  presence  of  an 
additional  nuisance  parameter  -  moment  of  signal  disappearance  7,  does  not  affect  the  final 
structure  of  the  test  procedures.  But  it  certainly  does  affect  the  detection  performance,  since 
the  accumulated  SNR  becomes  less  as  7  decreases. 

Another  possible  modification  of  the  algorithm  is  to  use  the  statistic  Uk  with  reflection  from  zero 
barrier, 

(6.12)  Uk  =  (Ak  +  Uk-1)+  ,  *>1,  #o  =  0, 

in  place  of  Uk  in  (6.10).  Obviously,  these  two  tests  have  the  same  performance  as  long  as  the 
threshold  is  positive.  Indeed,  it  is  easy  to  see  that  the  trajectories  of  the  statistics  Uk  and  Uk 
coincide  in  the  non-negative  half-plain. 

6.4.  Detection  Algorithm  2:  Fully  Sequential  Procedure.  We  now  change  the  problem 
set-up  and  consider  the  fully  sequential  procedure.  To  be  specific,  we  assume  that  k  =  n,  n- 1, . . . 
and  the  decision  about  target  presence  or  absence  is  made  at  each  moment.  This  case  is  probably 
the  most  relevant  to  the  IRST  surveillance  systems,  which  work  for  a  long  time  periodically  raising 
false  alarms. 

It  is  easy  to  see  that  the  generalized  LLR  L_oo,n  :=  Ln  is 


Ln  =  max(C/„,  l7»_i,-*-)> 


20 


CENTER  FOR  APPLIED  MATHEMATICAL  SCIENCES,  USC 


where  Uk  is  defined  in  (6.5).  Formally  the  sequential  algorithm  that  solves  the  problem  is 
identified  with  the  stopping  time 

ra  =  inf{n  :  Ln  >  o}. 

The  stopping  time  ra  can  be  rewritten  as 

(6.13)  r0  =  inf{n  :  Un  >  a}, 

which  is  exactly  the  CUSUM  procedure.  Indeed,  the  decision  on  target  absence  (dn  =  0),  which 
is  equivalent  to  observation  continuation  at  time  n,  is  made  when 

{Ln  Ln—1  ^  flj ■■  ■ }  {Un  ^  Un—i  ^  ■  " }  1 

Thus,  again  prior  uncertainty  with  regard  to  the  moment  of  disappearance  does  not  affect  the 
structure  of  the  detection  algorithm  (if  7  is  considered  as  a  nuisance  parameter)  but  it  does  affect 
its  performance. 

It  is  important  to  understand  that  the  algorithm  (6.13)  requires  additional  strategy  after  the 
decision  d  =  1  (target  appeared)  is  made.  Since  we  expect  multiple  signal  appearances,  one  may 
not  simply  stop  observations  after  this  decision  is  made.  One  way  to  handle  this  problem  is  to 
set  Un  =  0  when  ra  =  n  and  to  start  all  over  again  (immediate  renewal).  Then  the  algorithm 
is  immediately  ready  to  detect  the  next  target.  In  this  case  the  CUSUM  statistic  is  modified 
(compared  to  (6.5))  as  follows 

(6.14)  Un  =  An  ~b  Un_ill{o<{7n_i<a}> 

where  ll{y}  is  the  indicator  of  the  set  Y.  This  modification  has,  however,  one  drawback:  it  may 
(and  usually  does)  lead  to  multiple  detections  of  the  same  signal  and  hence  requires  additional 
identification  logic. 

Another  reasonable  way  is  to  use  the  detection  rule  (6.13)  with  pure  CUSUM  (6.5)  supplemented 
by  the  following  logic  for  new  target  detection.  Let  va  be  the  first  time  such  that  Un  goes  below 
the  threshold  a  after  exceeding.  To  be  precise, 

va  =  inf{n  >  ra:Un  <  o}. 

Then  for  ra  <  k  <  va  —  1,  one  confirms  the  presence  of  the  same  target,  while  if  for  some  m  >  1, 
UVa+m  >  a,  then  we  make  the  decision  that  a  new  target  track  appears  at  time  va  +  m.  In  other 
words,  if  the  i— th  target  was  detected  at  time  t, ft  and  confirmed  till  v®  - 1,  then  the  (i  +  l)-th 
target  is  detected  at 

ril+1)  =  inf {n  >  i/«  :Un>a} 
and  its  track  is  confirmed  till  the  moment  i4*+1)  -  1,  where 

ua+1^  ~  inf{™  >  Ta+V>  -Un  <  a}. 

We  also  note  that  one  may  replace  the  statistic  (6.14)  by  the  statistic 

(6.15)  Un  =  (An  +  Un-l)ll{0<An+t/n_i«i}  • 

This  replacement  does  not  affect  performance. 


DETECTION  ALGORITHMS  AND  TBD  ARCHITECTURE  FOR  IRST 


21 


Applying  a  similar  argument  to  the  quasi-Bayesian  statistic  Gn,  we  obtain  that  the  stopping  rule 

tb  =  inf {n  :  Gn  >  B} 


can  be  rewritten  as 

(6.16)  rB  =  inf{n  :  Rn  >  B}. 

In  the  case  of  immediate  renewal  the  average  likelihood  ratio  Rn  is  computed  recursively  according 
to  the  formula 

(6.17)  Rn  =  exp(A„)  (l  +  Rn-i l^-i <b})  • 

The  same  logic  as  above  can  be  used  to  detect  not  only  the  fact  of  target  appearance  but  also  its 
disappearance.  However,  in  the  case  of  multiple  targets  it  is  better  to  consider  a  more  complex 
set  of  hypotheses  that  would  include  at  least  three  alternatives:  the  last  signal  did  not  appear, 
appeared  but  is  still  present,  appeared  and  disappeared.  This  case  is  considered  in  the  next 
section. 

We  stress  once  again  that  the  statistic  Rn  may  be  interpreted  as  the  average  likelihood  ratio  over 
the  uniform  distribution  of  the  A  on  the  interval  [1,  n],  while  the  exp(U„)  is  the  maximum  of  the 
same  likelihood  ratio  over  A  €  [1,  n]. 

The  method  (6.16)  has  an  advantage  over  the  generalized  CUSUM  algorithm  -  the  statistic 
Rn  =  Rn  -  n  is  a  zero-mean  martingale  with  respect  to  P0  regardless  of  the  i.i.d.  assumption 
on  the  observations.  As  a  result,  if  one  sets  B  =  1/Er,  then  the  frequency  of  false  alarms  is 
upper  bounded  by  the  prespecified  value  Fr  (see  Section  6.7.1  for  details).  In  other  words,  the 
threshold  B  is  easily  estimated  for  a  large  class  of  models.  There  is  no  similar  result  for  the 
CUSUM  statistic  Un- 


6.5.  Algorithm  3:  Joint  Detection  of  Track  Appearance  and  Disappearance.  Consider 
an  extended  decision  making  process  assuming  that  after  each  decision  on  target  disappearance 
the  previous  data  are  discarded.  Specifically,  we  accept  the  following  logic: 

1.  If  the  decision  that  the  i-th  target  disappeared  is  made  (d,  =  DA),  then  one  of  the  two 
decisions  may  be  made  -  the  (i  +  l)-th  (new)  “target  did  not  appear”  (di+1  =  NA  =  0)  or 
the  (i  +  1)— th  “target  appeared  and  is  still  present”  (di+i  =  A&P). 

2.  If  the  decision  di+ 1  =  1  that  the  (i  +  1)— th  target  appeared  is  made,  then  one  of  the 
two  decisions  may  be  made  -  “target  is  still  present”  (di+1  =  P)  or  “target  disappeared” 
(di+ 1  =  DA). 

The  statistics  Ln  and  Gn  are  then  the  likelihood  (relative  to  H0  -  target  is  absent  at  all)  of  the 
composite  (and  combined)  hypothesis  H  that  includes  two  sub-alternatives: 

“ H  :  target  appeared  and  is  still  present  (Ha&p)  +  target  appeared  and  disappeared  (Fa&d)”- 
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Clearly,  at  time  moment  n  the  likelihood  of  the  subalternative  T/a&p  is  determined  by  the  statistic 
Un,  while  the  likelihood  of  the  subalternative  #a&d  by  Ln_x  =  maxfc<n_i  Uk.  Thus  the  structure 
of  the  decision  making  algorithm  may  be  as  follows. 

After  the  decision  d*  =  DA  (the  i-th  target  disappeared)  is  made,  the  statistic  Un  is  formed 
with  the  null  initial  condition.  This  statistic  is  compared  (at  each  step)  with  the  threshold  a.  If 
Un  <  a,  the  observation  is  continued.  Otherwise  ( Un  >  a),  the  decision  di+ 1  =  1  (the  (i  + 1)— th 
target  appeared)  is  made,  the  statistic  Un  is  computed  further  and  also  the  algorithm  starts  to 
compute  the  statistic  Ln  with  the  initial  condition  L  «+u  =  Uu+d,  where  riI+1^  is  the  moment 
of  detection  of  the  ( i  +  l)-th  target.  The  difference  A„  =  L„_i  -  Un  is  compared  with  another 
threshold  b  at  each  step  n  =  ri*+1)  +  1,  ril+1)  +  2, . . . .  If  An  <  b,  then  the  decision  di+i  =  P 
(still  present)  is  made  and  the  next  time  step  is  analyzed.  If  for  some  n  =  v^+l\  An  >  b,  then 
the  decision  di+i  -  DA  (target  disappeared)  is  made,  the  statistic  Ln  is  not  computed,  but  the 
statistic  Un  is  computed  for  n  =  v^+l)  +  1,  +  2, . . .  with  the  initial  condition  17  (<+i)  =  0. 

b 

Then  the  whole  cycle  is  repeated. 

Thus,  the  track  appearance  of  the  (i  +  1) — th  target  is  detected  at  the  moment 

t^+1)  =  inf{n  >  :Un>a} 
and  its  disappearance  is  detected  at 

i 'f+1)  =  inf{n  >  rW  :  A n>b}  . 

The  value  of  the  threshold  b  may  be  computed  based  on  the  trade-off  between  the  error  proba¬ 
bilities  due  to  too  early  decision  about  target  disappearance  and  the  delay  of  this  decision  when 
the  target  really  disappears. 


Figure  4.  Block  diagram  of  Algorithm  3 
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Figure  5.  Detection  of  track  appearance/disappearance  by  Algorithm  3 


The  similar  decision  making  algorithm  may  be  built  upon  the  statistics  Rn  and  Gn.  Then  the 
appearance  of  the  {i  +  l)-th  target  is  detected  at  the  moment 

fg+1^  =  inf{n  >  v®  :  log  Rn  >  log/?} 

while  its  disappearance  is  detected  at  time 

p(t+i)  _  jnf|n  >  fW  :  An  >  c}  , 

where  An  =  log Gn-\  -  log/4  with  the  initial  condition  A =  0  and  where  c  is  a  threshold. 

B 

The  block  diagram  of  the  algorithm  and  the  typical  behavior  of  the  decision  statistics  are  shown 
in  Figure  4  and  Figure  5,  respectively. 

6.6.  Adaptive  Detection  Algorithms.  Recall  that  in  general  the  position  of  the  target  is 
unknown  and  the  density  pi{Yk  \  9k)  depends  on  the  parameter  9k  that  characterizes  the  position. 
Then  instead  of  unknown  9k  one  may  use  the  estimate  of  the  position.  This  estimate  can  be 
obtained  by  applying  the  TbD  procedure  (based  on  ONF).  Thus  the  development  of  adaptive 
versions  of  the  above  detection  algorithms  is  needed. 

The  statistics  U„  =  ...,9n)  and  Rn  =  Rn(9u  ...,9n)  are  the  functions  of  the  sequence  of 

theta’s  till  time  n.  Thus  the  developed  procedures  can  not  be  applied  directly.  If  the  9n  can  be 
estimated  (this  is  the  case  in  our  system),  then  the  natural  solution  is  to  use  the  statistics  Un  = 
[4(0i, . . . ,  0„)  and  Rn  =  Rn(9\, . . . ,  0„),  which  are  computed  based  on  the  recursive  formulas 

Un  =  K(0n)  +  Uti,  /4  =  eAn(4)(l  +  £n-i). 

Here  9k  =  9k(Yf)  is  an  estimate  of  9k  based  on  the  previous  k  observations. 

However,  it  is  very  difficult  to  evaluate  the  performance  of  the  corresponding  adaptive  algorithms. 
The  reason  is  that  Un  is  not  a  CUSUM  and  Rn  is  not  an  average  likelihood  ratio  anymore.  In 
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fact,  the  exp{A„(0„)}  is  not  the  partial  likelihood  ratio,  since  pi(Yk  |  6k(Y ?))  is  not  a  probability 
density.  Particularly,  the  most  important  question  of  how  to  choose  the  thresholds  remains  open. 

To  avoid  this  difficulty  we  propose  the  following  trick.  At  stage  k  instead  of  using  the  estimate 
6k  we  propose  to  use  the  one  stage  delayed  estimate  0fc_i(Yifc_1).  In  other  words,  instead  of  Un 
and  Rn  we  will  use  the  statistics  U*  and  R*n,  which  satisfy  the  recursions 

(6.18)  U*  =  A*  +  [U*_x]+ ,  n  =  2,3,...,  U{  =  0; 

(6.19)  R*n  —  eAn(l  +  Rn~i)i  n  =  2,3,...,  R\  —  0, 

with  A*  =  An(0n-i).  The  corresponding  adaptive  target  detection  procedures  have  the  form 

(6.20)  r*  =  inf{n  >  2  :  U*  >  a}, 

(6.21)  f*B  =  inf{n  >2  :  R*n>B}. 

Since  pi(Yn  |  (9„_i)  is  the  probability  density  for  any  F"-1 -measurable  estimate  of  0n,  A*  is  the 
log-likelihood  ratio.  As  a  result,  the  statistics  (6.18),  (6.19)  preserve  most  of  the  nice  properties 
of  the  former  statistics  Un  and  R^.  In  particular,  R*n-n  is  a  P0-martingale  with  mean  zero. 
This  fact  allows  us  to  upper  bound  the  frequency  of  false  alarms  in  the  adaptive  algorithms 
regardless  of  specific  pre-change  and  post-change  distributions  (see  Section  6.7.1). 

The  block-diagram  of  a  typical  adaptive  track  appearance/disappearance  detection  algorithm 
that  uses  the  estimates  of  target  location  from  the  ONF-based  TbD  scheme  is  shown  in  Figure  6. 


Figure  6.  Block  diagram  of  a  typical  adaptive  detection  algorithm 


6.7.  Choice  of  Thresholds  and  Performance  Evaluation.  In  this  section  we  give  an  ar¬ 
gument  and  some  ideas  that  can  be  used  for  performance  evaluation  of  the  proposed  detection 
algorithms.  We  focus  on  the  first  two  algorithms.  The  analysis  of  the  third  algorithm  is  more 
difficult  and  we  leave  this  problem  for  the  future. 

In  both  algorithms  (Algorithm  1  and  Algorithm  2)  it  is  desirable  to  fix  the  probability  of  false 
alarm  P{&  in  a  given  fixed  interval  of  length  T.  We  note  that  while  in  Algorithm  1  this  is  the 
length  of  the  sliding  window,  in  Algorithm  2  this  is  the  “artificial” ,  auxiliary  parameter  that  is 
defined  based  on  tactical  conditions4.  For  the  fully  sequential  algorithms  (Algorithm  2)  it  is  also 


4If  the  dense  flow  of  targets  from  a  particular  direction  is  expected,  T  should  be  chosen  much  less  than  the 
typical  average  duration  of  a  signal.  Otherwise  T  should  be  comparable  with  the  average  time  of  target  presence. 
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reasonable  to  fix  the  frequency  of  false  alarms,  Fr  =  I/Eqt,  at  the  specified  level  Fr.  In  the  next 
subsection  we  show  that  it  may  be  done  in  the  general,  not  necessarily  i.i.d.  case. 


6.7.1.  Upper  Bounds  for  the  Frequency  of  False  Alarms.  Consider  the  general  case 
where  observations  can  be  dependent  and/or  non-stationary.  To  be  specific,  assume  that  under 
Ho,  p(Yi)  =  poOT)>  while  if  the  target  appears  at  the  moment  A  and  is  still  present  at  time  n 
the  model  is 

PA(y1n)=Po(nA_1)xpi(>TinA"1), 


where  p0(-)  is  the  joint  density  of  the  noise  and  pi(-)  describes  probabilistic  properties  of  the 
mixture  of  signal  and  noise.  By  ^(Y")  =  denote  the  likelihood  ratio  of  the  hypotheses 

H\  and  H0.  Obviously, 


9l 


9tl 


9n 


where 


9kx 


1 if'1) 

Po <Y}  |  Y?~l) 


Hence  for  the  CUSUM  and  the  average  likelihood  ratio  statistics  we  have 

(6.22)  Un  =:  max  log  =  max(0,  max :  log#”-1)  +  log#"  =  logp"  +  U+_x,  U0  =  0, 

\<n  A<n~l 

n  n— 1 

(6.23)  Rn  ='■  =  9n  X/5"-1  =  +  ■^n^’  ^0  =  0. 

A=1  A=1 


Particularly,  in  the  i.i.d.  case  log  #"(7”)  =  A n(Yn)  and  the  relations  (6.22)  and  (6.23)  coincide 
with  (6.2)  and  (6.3),  respectively.  In  the  latter  case  both  statistics  are  also  Markov  processes, 
which  helps  to  evaluate  performance  of  the  detection  algorithms.  But  in  general  the  statistics 
are  non-Markov  and  the  analysis  is  more  difficult. 


The  most  important  question  is  how  to  choose  the  thresholds  in  order  to  guarantee  a  specified 
false  alarm  rate.  To  answer  this  question  we  first  observe  that  -E0(#”  |  Y"-1)  =  1  for  any  n  >  1 
and  any  model.  As  a  result,  it  is  easily  seen  that  for  every  n  >  1  the  statistic  Rn  -  n  is  a 
jPo — martingale  with  mean  zero, 

E0(Rn  -  n  I  Yi"-1)  =  Rn—i  -  (n  -  1),  E0Rn  =  n, 


and  this  is  valid  for  an  arbitrary  model. 

Then,  by  the  optional  stopping  theorem,  for  any  Markov  time  r,  E0Rr  =  E0r,  while  by  the 
definition  of  the  stopping  time  tb,  E0RfB  >  B.  This  implies 

(6.24)  E0tb  >  B, 

which  holds  for  any  model  of  signal  and  noise.  Therefore  if  we  set 

(6.25)  B  =  1/PF, 

then  the  frequency  of  false  alarms  will  be  upper  bounded  by  Fr. 


26 


CENTER  FOR  APPLIED  MATHEMATICAL  SCIENCES,  USC 


It  turns  out  that  the  martingale  property  of  the  statistic  Rn  can  be  used  to  obtain  the  upper 
bound  for  the  frequency  of  false  alarms  in  the  CUSUM  procedure.  Indeed,  for  any  n  >  1 

n 

exp(?7n)  =  ma xjJ  <Rn  =  ^9x 

— n  A=1 

and  hence  on  the  event  {ra  <  oo},  exp(t/Ta)  <  Rra.  Thus,  if  E0ra  <  oo,  then 

ea  <  E0exp(UTa)  <  E0Rra  =  E0Ta  , 

where  the  left  inequality  follows  from  the  definition  of  the  stopping  time  ra  and  the  right  equality 
from  the  martingale  property  of  Rn.  This  shows  that  for  any  model  (not  necessarily  i.i.d.)  that 
guarantees  finiteness  of  E0ra  we  have  the  inequality 

(6.26)  E0ra  >  e“  . 

Hence  if 

(6.27)  a  =  log(l/Fr), 

then  the  frequency  of  false  alarms  in  the  CUSUM-type  procedure  is  upper  bounded  by  Fr. 

Finally,  we  note  that  (6.24)  and  (6.26)  also  hold  for  the  adaptive  sequential  algorithms  (6.20), 
(6.21),  since  these  algorithms  are  particular  cases  of  (6.22),  (6.23)  with  g "  =  exp{A*}. 

6.7.2.  The  Case  of  i.i.d.  Observations.  We  first  observe  that  in  the  i.i.d.  case  the  statistic 
Un  is  a  Markov  process  with  the  transition  probability  densities  (under  Hi,  i  =  0,1), 

(6.28)  pi(un\un_})  =  fi(un  -  g(un_ i)), 

where  fi(y )  is  the  probability  density  of  the  LLR  A k  =  Uk~  p(C4-i)  under  Hf,  g(y)  =  y+  for 
Algorithm  1  (fixed  sliding  window)  and  g(y)  =  yl{o<y<a}  f°r  Algorithm  2  with  immediate  renewal 
(see  (6.14)).  So  is  the  statistic  rn  =  log  Rn  with  the  transition  probability  density 

Pi(un  |  rn—i  —  Un—\)  fi^Un  lOg (l  "h  C  l{un_l<6})^ 

for  the  quasi-Bayesian  procedure  (see  (6.17)).  The  analysis  of  these  algorithms  may  be  done  by 
using  the  renewal  theory  (see,  e.g.,  [20,  46,  59]). 

In  what  follows  we  consider  only  the  CUSUM-type  procedure.  For  the  quasi-Bayesian  algorithm 
all  formulas  are  similar.  The  false  alarm  probability  is 

(6.29)  Pfa(T,  a)  =  Po{ra  <  T)  =  1  —  Pq{U\  <  a, . . . ,  Ut  <  a), 
while  the  target  missing  probability  is 

(6.30)  Pmis(T,7>  A)  =  P\n{ra  >  T)  =  P\n{Ui  <  a, . . .  ,UT  <  a), 
where  the  statistic  Un  obeys  the  recursion 

(6.31)  Un  =  An  +  g(Un- 1),  n  =  l,...,T. 

The  threshold  is  chosen  from  the  equation 

Pfa(T,  a)  =  a0, 
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where  ao  is  a  specified  false  alarm  rate.  Therefore  to  choose  the  thresholds  and  to  evaluate  the 
probabilities  of  correct  detection  one  needs  to  find  the  probability  of  entry  in  the  set  [a,  oo)  by 
the  Markov  process  Un  during  the  time  T.  This  process  starts  from  0  for  Algorithm  1  and  from 
some  random  point  in  (— oo,  a)  for  Algorithm  2. 


In  general,  let  {Xn,0  <  n  <  T},  X0  =  x  (x  is  fixed)  be  a  Markov  process  with  the  transition 
density  f(y,  z)  (Xn_i  =  z  -*•  Xn  =  y).  Define  the  system  of  functions  pk(y,  x)  by  the  recursion 


(6.32) 


Pk+i{y,x) 


Pk(y,z)f(z,x) 


dz, 


k  >  1 , 


Pi(y,x)=F(y,x) 


/y 

f(z,  x )  dz  . 

-oo 


The  function  pk(y,x)  is  nothing  but  the  conditional  probability  of  the  event  {Xk  <  y,Xk_ i  < 
a,...,Xi<a}  given  X0  =  x.  Then 


(6.33) 


P(Xi  <  a, . . . ,  XT  <  a  |  XQ  —  x)  =  pT(a,  x ) . 


6.7.3.  Performance  of  Algorithm  1.  Similar  to  (6.32)  for  i  =  0,1,  define  the  functions 
Pk\y,x)  by 


Pkhiy,  ®)  =  /  z)fi(*  -  x+)  dz>  k  ^ 1  > 

(6.34) 

Pi\y,x)  =  Fi(y-x+) , 


where  fi(y  -  x+)  is  the  transition  density  function  of  the  statistic  Un  (under  Hi)  that  satisfies 
(6.31)  with  g(Un-i)  =  and  the  null  initial  condition,  U0  =  0.  The  probabilities  of  errors  are 
then  computed  as  follows 

(6.35)  Pfa(T,  a)  =  1  -  pP  (a,  0) ; 

(6.36)  Pmis(T,A,7,a)=P?)(«,0),  if  A  <  n  -  T  +  1,  7  >  n . 


IfA>n-T+l  or/and  7  <  n,  then  the  situation  is  more  delicate.  It  may  be  shown  that  in  this 
case  the  probability  of  target  missing 


(6.37) 


■PmisCf)  A,7,  a)  =  f  f  pi°l7(a,  OpP-x+M^^Px-Mv,  °) 

J  —00  J  —00 


A  >  n-T  +  1,  7  <n 


To  obtain  the  mis-detection  probability  for{A<n  —  T+l,7<n}  and  {A>n  —  T+l,7>n}, 
one  has  to  put 

P\-i(dVi  0)  =  drl  ^or  {A  <  n  -  T  +  1, 7  <  n}  , 

Pnly(a,  £)  =  1  for  {A  >  n  -  T  +  1, 7  >  n} 

in  (6.37)  (S(x)  being  a  delta-function).  Particularly,  this  yields  the  relation  (6.36). 
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It  is  worth  mentioning  that  the  probabilities  p%\a,  0)  satisfy  the  integral  equations 

PfcV,  0)  =  [  p(Z)_l{a,z)fi(z)dz,  k>  1 , 

(6.38)  ■'-oo 

Po\y,x)  =  1  > 

which  are  simpler  than  (6.34). 

The  formulas  (6.35)-(6.37)  determine  the  performance  in  terms  of  the  probabilities  of  false  alarm 
Pfe  and  target  missing  Pmis  at  each  current  time  n  within  the  window  of  the  length  T.  In 
particular,  the  threshold  a  is  found  from  the  equation  (see  (6.29)) 

1  -  a0  =  pP  (a,  0) 

where  the  function  p?\a,  0)  satisfies  the  recursive  relation  (integral  equation)  (6.38).  Also,  the 
probability  of  target  missing  in  the  interval  of  its  presence,  [A,  7],  can  be  approximated  by5 

P mis(7  —  A  +  1,  d)  ~  ifli  0)  • 


The  following  elementary  estimates  are  useful  for  initial  choice  of  threshold  and  target  missing 
probability  evaluation: 

(6.39)  Pfa(T,  o)  >  P0(ST  >  a),  Pmis( 7  -  A  + 1,  a)  <  1  -  Pa,7(S7  -  SX-i  >  a) 
where  Sn  =  YYk=i  A*. 


Along  with  the  probabilities  of  errors  it  is  interesting  to  evaluate  the  mean  time  delay  in  target 
detection  PA]7{ra(T,  n)  —  A|ra(T,  n)  >  A}  (in  order  to  estimate  benefits  of  the  sequential  scheme 
(6.10)  compared  to  the  direct  non-sequential  test  (6.8)).  Here  E\tl  is  the  symbol  of  expectation 
for  the  given  values  of  A  and  7.  By 


'piivY 

Mv). 


Pi(y)dy 


denote  the  Kullback-Leibler  information  number.  In  case  when  T  a/ 1,  7  —  A  >  a/ 1,  and  a  is 

sufficiently  large  (cn0  is  small)  the  mean  detection  delay  may  be  approximated  by 


E\ )7{r0(T,n)  -  A| rft(T,n)  >  A}  «  a/ 1. 


6.7.4.  Performance  of  Algorithm  2.  For  Algorithm  2,  reasonable  performance  may  be  ex¬ 
pressed  in  terms  of  the  frequency  of  false  alarms  Fr  =  1/E0ra  ( E0ra  is  the  mean  time  between 
false  alarms)  and  the  detection  delay  E\^{ra  -  A|r0  >  A}.  However,  such  characteristics  as  the 
probabilities  of  false  alarm  and  target  detection  in  the  fixed  interval  of  the  length  T  are  also  of 
interest. 

For  Algorithm  2  with  immediate  renewal  we  have 

(6.40)  Fr  <  exp(-a),  ExnM  ~  ^lTa  >  A}  »  a/I. 

5 This  formula  is  approximate,  since  U\-i  is  a  random  variable  belonging  to  the  interval  (-00,  a).  For  A  =  1  it 
is  exact. 
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The  first  formula  in  (6.40)  follows  from  (6.26)  and  the  second  formula  is  valid  when  7  -  A  >  a/I 
and  a  is  sufficiently  large  (Fr  is  small).  According  to  the  first  inequality  in  (6.40)  the  threshold 
may  be  chosen  as  a  =  log(l /Fr)  to  guarantee  the  frequency  of  false  alarms  Fr  <  Fr,  where  Fr  is 
a  specified  level. 

However,  it  is  even  more  reasonable  to  choose  the  threshold  based  on  the  specified  probability 
of  false  alarm  Pfa(T,  a)  within  the  time  interval  of  the  length  T  at  each  current  point  in  time  n. 
Then  the  mean  detection  delay  may  be  approximately  evaluated  by  the  second  equality  in  (6.40) 
and  the  algorithm  is  optimal  -  it  minimizes  the  mean  detection  delay  among  all  algorithms  with 
this  kind  of  constraints  (at  least  for  sufficiently  large  a). 


To  obtain  the  probabilities  of  errors  in  the  fixed  interval  T  we  note  that  the  statistic  Un  is  again 
the  Markov  process  which  obeys  the  recursion  (6.31)  with  the  random  initial  condition  U0  =  u 
and  with  g(Un- 1)  =  tIn-il{o<t/n_i<a}-  The  random  variable  u  has  the  distribution  G0(y),  which 
is  nothing  but  a  stationary  distribution  of  the  Markov  process  {Un}  under  the  hypothesis  H0 
(target  is  absent).  This  distribution  satisfies  the  integral  equation 


G0(x)  —  J'  Fq(x  sH{o<s<a})  d(?o(s) , 


where  F0(y)  is  the  distribution  function  of  An  under  H0.  Thus,  the  false  alarm  probability  is 
found  as 

Pfa(T,  a)  =  (  Pfa  (T,  a,  x)  dG0  (x) , 

J  —OO 

where  Pfa(T,  a,  x)  is  the  probability  (6.29)  in  case  when  U0  —  x,  x  <  a.  This  last  probability  is 
computed  according  to  (6.33)  with  pr(a,  x)  =  pP(a,  x),  which  is  defined  by  the  recursive  relation 

/a 

Pfc— 1(®>  ^O/o  (2  —  ®l{o<*<a})  dz,  k  >  1 ,  p[  \a,x)  =  1  . 

-00 


As  we  mentioned  above,  the  analysis  of  the  proposed  algorithms  may  be  done  based  on  the 
renewal  theory.  This  theory  allows  for  more  accurate  performance  evaluation  as  compared  to 
the  above  formulas  (for  the  mean  detection  delay  and  frequency  of  false  alarms).  The  detailed 
analysis  of  the  proposed  algorithms  and  their  tuning  will  be  done  in  the  near  future  (see  also 
Section  6.8,  relationships  (6.41), (6.42)  for  more  accurate  estimates  in  a  Gaussian  case). 


6.8.  TbD  Model  with  Spatio-Temporal  Matched  Filter.  Let  us  apply  the  above  results 
to  the  case  where  the  target  moves  along  a  deterministic  trajectory.  Then  the  spatio-temporal 
matched  filter  is  the  optimal  tool  for  TbD.  Since  target’s  velocity  and  direction  of  motion  are 
unknown,  the  bank  of  filters  should  be  used  and  the  detection  is  done  in  each  channel  indepen¬ 
dently.  Assume  also  that  after  clutter  suppression  the  residual  clutter  (plus  noise)  samples  are 
i.i.d.  and  Gaussian,  Af(0,(j2)  (see  the  model  (5.9)).  Then  the  LLR  A*  is 

Afc  =  “ 2  S  ^  -  Xii  ~(k~  1)vA)Za(k)  ~(k~  1>A)2  . 

i,j  id 

In  the  above  formula  we  implicitly  used  the  assumption  that  the  target  moves  with  the  constant 
given  speed  v  and  in  a  given  direction  (i.e.,  in  the  particular  “velocity-angle”  channel).  This 
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formula  is  readily  generalized  for  arbitrary  (but  known)  dynamics.  Note  that  here  we  use  the 
notation  Yk  =  {Zij(k)}. 

According  to  (6.31)  the  statistic  Un(x)  is  calculated  as 

Un(x)  =  An(x)  +  g(Un-i(x)),  n  =  1, 2, . . . 

Then  the  above  algorithms  are  applied  in  each  channel  independently  (here  we  mean  the  channels 
that  are  related  to  different  values  of  the  initial  position  x). 


If  the  signal  is  exactly  located  in  the  given  channel,  then  the  LLRs  {A* (a)}  are  i.i.d.  Gaussian 
random  variables  with  parameters 

EiAk  =  -Eq  Ak  =  \p,  Var0Afc  =  VariAfe  =  \p, 
where  p  =  ^  Yhi,j  h{xij)2  *s  SNR. 

Performance  of  the  developed  algorithms  is  evaluated  by  formulas  obtained  above  in  Section  6.7 
with  the  densities 

My)  =  If  (fTf) '  m  =  ff  ' 

where  <p(y)  is  a  standard  normal  density  function, 

^  =  v^exp{~^} ' 

In  particular,  for  Algorithm  1  (sliding  window),  the  simplest  estimates  for  probabilities  of  errors 
are  given  by  (see  (6.39)) 

2a  -  p( 7  -  A  +  1) 


Pfa(T,a)>l  —  Pmis(7-A  +  l,a)<$ 

where  $(y)  is  a  standard  normal  distribution  function. 


y/piri-  A  +  l) 


Also,  for  the  sequential  Algorithm  2  with  immediate  renewal  for  sufficiently  large  a  and  7  >  2 a/ p 
the  following  approximate  estimates  hold: 


(6.41) 

where 


(6.42) 


Ft  <  0e-‘  , 

-  ^|t«>A}as^(a  +  >t  +  C), 

^  =?exp|-2E£a*(-lv^)|i 
=  El  [min{0,minn>iX)Li  Afc}]  e 


These  formulas  were  obtained  by  using  the  nonlinear  renewal  theory  [59].  They  improve  the 
upper  bounds  for  the  frequency  of  false  alarms  and  the  estimates  for  the  mean  detection  delay 
presented  above  in  the  general  case  (see  (6.26),  (6.27),  (6.40)). 
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Instead  of  making  decisions  on  target  presence  in  each  channel  we  may  apply  the  algorithm, 
which  is  based  on  the  maximum  likelihood  statistic, 

Mn  =  max  C/n(a:),  n  =  1,2,... 

X 

This  statistic,  however,  does  not  allow  for  a  convenient  recursive  computation.  To  simplify  the 
algorithm  we  use  the  “trick”,  which  was  already  used  in  Section  6.6  to  build  adaptive  detection 
algorithms  in  general  case.  By  x*n  —  x*n(Zn)  denote  the  maximum  likelihood  estimate  of  the 
target  location,  i.e. 

n 

x*n  =  argmaxV'  Afe(a;). 

X  Z J 

k= 1 

The  decision  statistic  is  defined  by  the  recursive  formula 

^  =  An«-1)+^:-i),  n  =  2,3, . . .  ,  U{=  0. 

The  “stopping  rule”  is  then  defined  by 

t*  =  inf {n  >2  :U*>a}. 

Thus,  we  use  the  one-stage  delayed  estimates  of  the  signal  location  in  the  partial  LLR  An.  The 
threshold  is  chosen  according  to  the  relation  (6.27),  which  allows  us  to  upper  bound  the  false 
alarm  rate.  This  algorithm  is  especially  attractive  when  only  a  single  target  is  expected.  In  the 
multi-target  situation  the  previous  algorithm  (decisions  in  independent  channels)  is  perhaps  the 
best  one  can  do. 

6.9.  TbD  Model  with  Optimal  Nonlinear  Filter.  In  the  context  of  TbD  with  optimal 
nonlinear  filtering  (see  Section  5)  the  estimate  of  the  partial  LLR  An  obtained  on  the  basis  of 
previous  n  —  1  observations  is  defined  by 

(6.43)  An  =  'y  ^  hn-i(xjj) Zjj (k)  —  —  ^  ^  hn—\{xjj)  , 

i,j  hi 

where  hk{xij)  =  h(Xk  -  Xij)  and  Xk  is  the  estimate  of  the  signal  location  based  on  the  k 
observations.  The  latter  estimate  is  formed  at  the  output  of  the  optimal  nonlinear  filter  that 
tracks  the  target  location.  We  will  consider  two  estimates  -  the  mean  estimate  and  the  maximum 
a  posteriori  estimate,  which  are  defined  below  in  (7.3)  and  (7.4)  (see  Section  7).  Thus,  the  above 
results  are  applied  by  using  the  LLR  estimate  (6.43). 

To  be  precise,  the  adaptive  CUSUM  statistic  is  defined  by  the  recursion 

(6.44)  U*  =  A*n  +  g{U*n_x),  n  =  2,3,...,  U*,  =  0 

and  the  “stopping  rule”  is 

(6.45)  t*  =  inf{n  >  2  :  [/*  >  a}. 

In  other  words,  again  we  have  used  the  same  trick  -  the  true  LLR  An  at  stage  n  is  replaced  with 
its  estimate  A*  =  An(Xn_i)  obtained  based  on  the  previous  n  —  1  observations.  Instead  of  Xn_i 
one  may  use  the  one  step  predicted  (extrapolated)  estimate  X*. 
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The  alternative  algorithm  (see  (6.19),  (6.21))  is  of  the  form 

(6.46)  t*b  =  inf{n  >  2  :  R*n  >  B}, 

where  the  estimate  of  the  average  likelihood  ratio  R*n  satisfies  the  recursion 

(6.47)  i?;  =  exp{A;}(l  +  J?;_1l{fl;_1<B}),  n  =  2,3,...,  R\  =  0. 

The  results  of  application  of  these  algorithms  will  be  discussed  in  Section  7.4. 


7.  Testing  of  the  Developed  Algorithms  for  IRST  Data.  Results  of  Experiments 

and  Simulation 

In  the  examples  considered  below  because  of  the  uncertainty  of  possible  behavior  of  a  non- 
cooperative  target,  the  alignment  of  successive  frames  is  extremely  difficult.  The  frame  alignment 
is  done  recursively  along  the  maximum  posterior  distribution  path  on  the  basis  of  the  optimal 
nonlinear  filtering.  Use  of  nonlinear  algorithms  is  necessitated  by  the  specifics  of  the  observation 
structure.  Indeed,  in  TbD  the  observation  function  (signal)  S(k,r)  is  highly  sharp  (essentially  a 
delta-function  in  small  domain  like  2  x  2  or  3  x  3  pixels).  In  addition,  IR  imaging  sensors  provide 
angle-only  measurements.  These  types  of  measurements,  especially  in  low  S(N+C)R  situations, 
practically  rule  out  application  of  extended  Kalman  and  similar  filters.  As  already  mentioned 
above,  application  of  multi-dimensional  matched  filters  and  banks  of  velocity  filters  is  possible 
only  for  constant  speed  targets  with  linear  trajectories,  which  is  not  the  case  in  the  examples 
considered.  At  the  same  time  the  proposed  ONF-based  algorithm  works  perfectly  well  even  for 
very  low  S(N+C)R  (up  to  -6dB  after  clutter  removal).  Particularly,  it  may  be  seen  from  the 
figures  presented  below  that  the  uncertainty  is  completely  overcome  after  processing  of  several 
frames  (in  these  examples  S(N+C)R  ranges  from  -ldB  to  -7dB  after  simple  pre-processing  and 
clutter  removal). 

The  developed  adaptive  detection  algorithms  that  use  the  estimates  of  target  location  constructed 
on  the  basis  of  ONF  TbD  work  also  perfectly  well.  Tracks  of  randomly  appearing  and  disap¬ 
pearing  targets  are  reliably  detected  up  to  -7dB  and  the  algorithms  allow  us  to  obtain  low  false 
alarm  rate  even  in  a  heavy  cluttered  background. 

7.1.  Clutter  Removal:  Real  IRST  Data.  The  proposed  nonparametric  approach  to  clut¬ 
ter  removal  (see  Section  4.1)  has  certain  advantages  over  conventional  methods.  This  fact  is 
illustrated  by  Figure  7  which  presents  the  results  of  kernel  smoothing  of  IR  data.  The  ‘local 
mean  removal’  procedure  used  for  clutter  rejection  in  a  series  of  papers  by  Reed  et  al.  (see  [39], 
[38])  is  equivalent  to  a  kernel  estimator  with  a  product  kernel  K  generated  by  uniform  kernels 
Ki(x)  =  K?,{x)  =  0.5I{N<i}.  It  may  be  seen  that  the  Fuller  kernel  (see  [26])  provides  superior 
smoothing  performance  compared  to  the  uniform  kernel.  The  reason  is  that  our  approach  relies 
on  the  trade-off  between  a  squared  bias  and  a  variance  of  estimators  while  the  uniform  kernel 
minimizes  only  the  asymptotic  variance  of  estimators  and  hence  induces  a  large  bias  term. 
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(a)  IR  data:  target  plus  clutter  (b)  Kernels  in  [-1,1] 
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(c)  Residuals  after  uniform  kernel  (d)  Residuals  after  Fuller  kernel 


Figure  7.  NAVY  IRSS  sensor  -  the  LAPTEX  field  test.  Original  IR  image  (a) 
and  the  results  of  clutter  removal  by  uniform  kernel  (c)  and  Fuller  kernel  (d) 

7.2.  Example  1:  TbD  of  a  Target  Based  on  IRST  Data.  The  first  example  presented 
is  a  naive  model  of  an  actively  maneuvering  air  target  which  flies  overhead  over  a  stationary 
observer  situated  on  the  ground.  The  complete  field  of  view  of  the  observer  is  modeled  by 
a  square  [0, 1]  x  [0, 1]  with  the  horizon  line  at  y  =  0.1.  The  target  initially  appears  close  to 
the  horizon  at  a  point  uniformly  distributed  inside  a  rectangle  [0.2, 0.8]  x  [0.1,0.15],  then  the 
horizontal  coordinate  evolves  as  a  Wiener  process  with  variance  a2t,  and  the  vertical  coordinate 
as  an  exponent  e^. 

Formally,  the  trajectory  model  is  described  by  a  stochastic  differential  equation6 
X0  ~  Unif(0.2, 0.8),  Y0  ~  Unif(0.1, 0.15); 

dXt  =  adWu  a  =  0.05,  dYt  =  (3Ytdt,  /?  =  0.6. 

The  target  was  assumed  to  emit  a  3  x  3  pixels  signature  which  was  observed  with  additive 
Gaussian  white  noise.  More  precisely,  let  <5  be  the  size  of  a  pixel  in  the  observed  image  in  both 


6Here  the  components  of  the  vector  r  are  denoted  by  (x,y).  Also  X  —  {X.  Y). 
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horizontal  and  vertical  directions.  The  target  shape  is  given  by 

h(x,  y)  =  l{|x|<3«/2}  l{|;/|<3<5/2} ,  (x,  y)  €  R2 , 

and  the  observation  at  time  moment  tk  =  kA,  A  =  0.05  with  noise  of  standard  deviation  a  and 
SNR  =  20  log  J  is 

Zij(k)  =  h{XkA  -  (i  -  1)5,  rfcA  -  (j  -  1)6)  +  a^j{k),  &(*)  ~  Norm(0, 1)  i.i.d. 

Examples  of  a  target  trajectory  and  IR  data  images  are  shown  in  Figure  8  and  Figure  9. 

Several  experimental  results  follow. 

•  For  moderate  noise  (a  =  1.232,  SNR=  -1.8dB)  the  simplest  and  fastest  Haar  basis  works 
quite  well  (see  Figure  10). 

•  For  more  intense  noise  (a  =  1.581,  SNR=  -4.0dB),  the  algorithm  utilizing  Haar  basis  loses 
track  occasionally  but  later  recovers  it  (see  Figure  11). 

•  For  borderline  cases  (a  =  1.679,  SNR=  -4.5dB),  the  Haar  basis  no  longer  provides  good 
tracking,  however  smoother  wavelet  bases  of  the  same  resolution  (such  as  Coiflet-1)  still 
achieve  good  tracking  accuracy  (see  Figures  12  and  13;  both  simulations  were  run  on  exactly 
the  same  observations,  and  the  difference  in  the  quality  of  tracking  is  quite  obvious). 

•  With  larger  noise  (cr  =  2.109,  SNR=  -6.5dB),  the  optimal  filter  fails  to  track  the  target, 
no  matter  what  basis  is  utilized. 


Trajectory 


Figure  8.  Typical  trajectory  of  overhead  target 
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Observation  at  t=0 
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Figure  9.  Sample  of  a  single  frame  of  observation 
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State  and  estimates,  coord=1 


time 

Confidence  intervals,  x 


State  and  estimates,  coord=2 


Estimate  errors 


Figure  10.  ONF,  Haar  64  x  64  basis,  a  =  1.232,  SNR=  -1.8dB 
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Figure  11.  ONF,  Haar  64  x  64  basis,  a  =  1.581,  SNR=  -4.0dB 
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Figure  12.  ONF,  Haar  64  x  64  basis,  a  =  1.679,  SNR=  -4.5dB 
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Figure  13.  ONF,  Coiflet-1  64  x  64  basis,  a  =  1.679,  SNR=  -4.5dE 
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State  and  estimates,  coord=1 
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Figure  14.  ONF,  Haar  64  x  64  basis,  a  =  2.109,  SNR=  -6.5dB 
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7.3.  Example  2:  TbD  of  a  Surface  Skimming  Missile.  Assume  that  a  stationary  obser¬ 
vation  platform  (camera)  is  at  a  fixed  height  h2  above  the  sea  level,  and  a  target  approaches 
with  constant  velocity  v  at  constant  height  hi  <  h2  above  the  sea  level.  Denote  the  Earth  radius 
R  =  63, 750, 000m.  Projection  center  (the  focal  point  of  the  camera)  is  situated  at  the  point 
with  coordinates  (0,0,  R  +  hi).  At  time  moment  t  =  0  the  target  appears  at  the  horizon  at 
angle  (3  clockwise  from  the  direction  of  y-axis  and  moves  at  angle  7  to  the  line  of  sight  from  the 
observer.  Because  of  rotational  invariance  after  projecting  onto  the  observation  cylinder,  assume 
that  P  =  0. 


The  three  dimensional  coordinate  system  used  next  has  the  origin  at  the  center  of  the  Earth, 
the  2 -axis  passing  vertically  through  the  observer,  the  y-axis  passing  horizontally  through  the 
observer  in  the  direction  of  the  point  where  the  target  becomes  first  visible  on  the  horizon,  and 
a;— axis  oriented  so  that  the  xyz  coordinate  system  is  a  Cartesian  right  hand  coordinate  system. 


xQ  =  0, 

y°  =  ^  ^  (^\J^Rh\  +  ~h\  +  yj 2Rh2  +  h^j  , 

R2  ~  +  h^2Rh*  +  h» 


R  +  h\  R  +  hi 


vt 


xt  =  (R  +  h2)  sin  7  sin  R  +  fi~, 
vt 


vt 


Vt  =  y0  cos 


Zt  =  Zo  cos 


R  +  h<L 
vt 

R  H-  h>2 


—  zo  cos  7  sm 


+  y0  cos  7  sin 


R  +  h<i 
vt 

R  +  ^-2 


After  projecting  onto  an  observation  cylinder  with  radius  r  (the  focal  length  of  the  camera),  the 
equation  of  the  motion  in  two  dimensional  coordinates  ( 6 ,  zp)  are  as  follows 

rzt 


Xt 

Qt  =  arctan  — , 

yt 


7v  _ 

zt  — 


VW+yl 


The  observed  signal  was  assumed  to  be  a  3  x  3  pixels  target  with  additive  Gaussian  white  noise. 
More  precisely,  let  <5i  be  the  size  of  a  pixel  in  the  horizontal  (coordinate  6 )  direction,  and  82  be 
the  pixel  size  in  the  vertical  (coordinate  zp)  direction.  The  target  shape  is  given  by 

(7.1)  h(x i,x2)  =  l{|*i|<3«i/2}l{|x2|<3fi2/2}>  (®i>  ^2)  e  M2, 

and  the  observation  at  time  moment  4  =  kA  with  noise  of  standard  deviation  a  and  SNR 
=  20  log  \  is 

=  @min  "b  (®  1)^1  > 

=  zpmin  +  (j  -  1)6,, 

(7.2) 

Zij  (k)  =  h(6kA  -  Qi ,  zpkA  -  z?)  +  ofa (k) , 

£ij(k)  ~  Norm(0, 1)  i.i.d. 


It  is  important  to  emphasize  that  the  Gaussian  model  for  the  residual  clutter  and  noise  is  used 
only  in  ONF  and  detection  algorithm  development.  The  algorithms  were  tested  with  the  use  of 
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the  real  IR  background  obtained  from  SPAWAR.  Particularly,  the  figures  in  Section  7.4  below 
contain  the  results  of  testing  for  these  real  cluttered  images  (i.e.,  {£ij(k)}k> x  in  (7.2)  are  given 
realistic  frames)  with  artificially  inserted  targets  whose  motion  and  shape  correspond  to  models 
(7.1),  (7.2). 


7.3.1.  TbD  of  Several  Targets.  Ability  to  detect  and  track  several  independent  targets  was 
implemented  by  running  the  corresponding  number  of  filters  in  parallel. 

The  exact  updating  formulas  are  given  below. 

Translation  operators.  The  operator  translates  the  argument  of  a  function  defined  on  the 
state  space  along  the  trajectory  of  the  model  #i  Let  us  denote  the  trajectory  of  the  model 
starting  at  the  point  (x,y)  at  time  moment  t  =  0  by  qi{t,x,y),  and  its  inverse  with  respect 
to  the  variable  t  (with  variables  x  and  y  fixed)  by  f(g,x,y).  Let  the  function  to  translate  be 
/  :  [0,  a)  x  [0;  b ]  x  {1}  U  [0,  a)  x  [-b,  b]  x  {-1}  U  {0}  x  {0}  x  {0}  -*  R.  Several  cases  have  to  be 
split  in  order  to  describe  the  translation  operator: 

C^a/)  Vi  1)  =  0)j  f°r  xg[0,o),  y  G  [0,  q  (x, 0,  A)), 

Pa/)  WA>x,2/),l)  =  for  x  e  [°>a)>  V  e  [0,gl(x,6,-A)), 

=f(x,y,l),  for  x  €  [0,o),  y  6  -A), 6), 

(2l/)(9*(A,x,t/),-l)  =/(*,y,-l),  for  are  [0,  a),  ye[0,b), 

(7l/)(0,0,0)  =  e-^/(0,0,0)  +  ^  J*  ’  ’  f(x,y,-l)dydx. 

Updating  formulas.  The  unnormalized  filtering  densities  plt(x,  y,  u)  corresponding  to  the  model 
are  recomputed  based  on  observations  Zk  =  {Zij(kA)}  with  time  step  A  according  to  formulas 

Po(0, 0, 0)  =  1, 

Plo(x,y,u)=  0,  for  (x,y,u)  ^  (0,0,0), 


pLa(0,0,0)  =  (TlAp\k_1)A)  (0,0,0), 

HkA(x,  y,  u,  Zk)  =  exp{^  ^  h(x  -xuy-  yj:  u)Ztj(kA)  -  &  J2  h(x  “  Xi>  V  ~  w)2}> 

plkA{x,y,u)  =  HkA(x,y,u,Zk)(TlAp\k_1)A)(x,y,u),  for  (x,y,u)  ±  (0,0,0). 


Next,  the  filtering  densities  are  normalized  as  follows: 


*fcA(*.V»«)  = 


A  A  =  PfcA(°>°,0)  +  J  JAa(*,VA) dxdy  +  J  J  p‘kA(x,y,-l)dxdy, 
plkA(x,y,u) 


n 


kA 
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and  two  estimates  of  the  position  of  the  l- th  target  are  computed.  Namely,  the  mean-square-error 
estimate: 

XlkAean  =  f  f  xnlkA(x,  y,  1  )dxdy+  f  f  xnlkA(x,  y,  -1)  dx  dy, 

(7-3)  .  r  r  r  r  , 

■yrijnean  _  j  j  yirlkA(x,y,l)dxdy  +  J  J  ynlkA(x,  y, -1)  dx  dy, 
and  the  maximum  o  posteriori  estimate: 

(7.4)  (xSax,  ylk a“)  =  ars  max  4a  (*,  y,  «)  • 


Both  estimates  can  be  used  in  different  situations.  The  first  estimate  is  theoretically  optimal. 
However,  it  usually  takes  longer  to  converge  to  the  real  position  of  the  target.  The  second 
estimate  “detects”  the  target  much  faster  but  is  more  prone  to  lose  track  due  to  noise. 

Precomputing.  Computationally,  the  most  expensive  part  of  the  algorithm  was  interpolating 
the  values  of  (TlAplt)  on  the  shifted  along  the  trajectory  (and  therefore  no  longer  uniform)  grid 
back  to  the  original  uniform  grid.  To  speed  up  the  computations,  off-line  part  of  the  algorithm 
computes  a  (very  sparse)  matrix  A  such  that  multiplying  the  vector  of  values  on  the  non-uniform 
grid  by  A  produces  the  vector  of  values  on  the  uniform  grid. 

Two  interpolation  algorithms  were  implemented:  nearest  neighbor  interpolation  (then  matrix 
A  has  only  zero  and  one  as  elements,  and  each  row  has  at  most  one  non-zero),  and  linear 
interpolation  using  known  values  on  a  triangle  surrounding  a  point  on  the  new  grid.  Linear 
interpolation  proved  to  be  much  more  stable  and  precise  in  simulations. 


7.3.2.  Simulation  results.  The  simulation  producing  results  included  below  used  two  targets 
and  the  noise  with  standard  deviation  a  =  1.4  which  corresponds  to  SNR  =  — 2.9dB. 

•  In  Figure  17  (t  =  Is),  the  target  is  present  but  not  yet  detected. 

•  In  Figure  18  ( t  =  6s),  a  single  sharp  peak  has  formed,  and  it  follows  the  target. 

•  In  Figure  19  (i  =  11s),  the  second  target  just  appears  at  the  horizon,  algorithm  doesn’t 
track  it  yet. 

•  In  Figure  20  ( t  =  16s),  both  targets  were  detected  and  are  being  tracked. 

Also  we  refer  to  enclosed  files  (multiframe  GIF  file  a6.gif  or  MS  video  file  a6.avi)  to  see 
the  incoming  observation  frames  (Navy  IRSS  sensor  -  the  LAPTEX  field  test)  with  very  dim 
synthetic  target  inserted.  The  second  pair  of  files  (multiframe  GIF  file  a5 .  gif  or  MS  video  file 
a5 .  avi)  shows  the  locations  of  two  targets  being  tracked  and  the  filtering  density  which  develops 
prominent  peaks  that  follow  both  targets. 
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Figure  15.  Surface  skimming  missile  trajectory  derivation 


Trajectories 


Figure  16.  Surface  skimming  missile  trajectories 
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x10-"  t=1.0s,  SNR=-6.7dB  Target:  (1. 45, -2.50e-004) 


Figure  17.  Surface  skimming  missile  tracking,  t  =  Is 


Figure  18.  Surface  skimming  missile  tracking,  t  =  6s 
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Figure  19.  Surface  skimming  missile  tracking,  t  =  11s 


Figure  20.  Surface  skimming  missile  tracking,  t  =  16s 
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7.4.  Appearance /Disappearance  Detection  of  a  Skimming  Missile.  In  this  section  we 
present  the  results  of  application  of  the  detection  algorithms  developed  in  Section  6.  The  model 
is  the  same  as  in  Section  7.3.  The  results  of  TbD  (the  estimates  (7.3)  and  (7.4)  of  the  signal 
location)  are  used  in  the  LLR  estimate  (6.43).  In  the  fully  sequential  CUSUM  algorithm,  the 
adaptive  CUSUM  statistic  U*  is  computed  recursively  according  to  (6.44)  and  then  is  compared 
to  the  threshold  a  (see  (6.45)).  In  the  sequential  (adaptive)  quasi-Bayesian  algorithm  the  estimate 
of  the  logarithm  of  average  likelihood  ratio  r*  =  log  R*n  is  computed  according  to  the  recursive 
formula  (6.47).  It  is  then  compared  to  the  threshold  log-B  (see  (6.46)). 

Figure  21,  Figure  22,  and  Figure  23  illustrate  the  behavior  of  the  adaptive  CUSUM  statistic 
for  different  S(N+C)R  (SNR  =  -0.25dB,  SNR  =  -3.62dB,  and  SNR  =  -6.6dB,  respectively). 
Theoretically  0  <  r*  -  U*  <  logn.  We  have  not  observed  any  substantial  difference  between 
these  statistics  in  our  experiments  -  the  difference  was  always  negligible.  So  we  omit  the  plots 
of  the  statistic  r*  and  show  only  the  behavior  of  the  adaptive  CUSUM.  It  turns  out  that  the 
mean  estimate  (7.3)  provides  better  results  than  the  maximum  posterior  estimate  (7.4)  and  we 
plot  only  the  CUSUM  statistic  that  uses  the  former  one.  It  is  seen  that  most  of  the  time  the 
statistic  U*  is  close  to  zero  or  negative  when  the  target  is  absent.  Sometimes,  however,  peaks 
arise.  These  peaks  may  be  identified  with  signal  presence.  But  they  are  typically  short  and 
may  be  easily  distinguished  from  the  peaks  due  to  target.  One  possible  method  to  discriminate 
between  false  alarms  and  true  target  is  to  make  the  decision  on  target  presence  if  there  are 
several  subsequent  exceedances  (say  3-4)  of  the  threshold.  Otherwise  the  decision  on  target 
absence  is  made  (i.e.  a  single  exceedance  is  identified  with  false  alarm).  When  the  target  appears 
the  statistic  rapidly  increases  while  when  it  disappears,  U*  decreases.  In  our  experiments  we 
observed  visible  difference  in  behavior  of  U*  when  the  signal  appears  compared  to  the  case  of 
its  absence  up  to  the  SNR  -6.6dB.  In  the  pictures  the  first  target  appears  at  time  n  =  1  and 
disappears  at  n  =  28.  The  second  target  appears  at  time  n  =  39  and  never  disappears. 

Figure  24  shows  the  results  of  detection  of  track  appearance  and  disappearance  by  Algorithm  3 
for  SNR  =  —  0.25dB  and  SNR  =  — 6.6dB.  The  detection  of  tracks  occurs  when  the  adaptive 
CUSUM  U*  exceeds  the  threshold  a  (the  upper  one)  and  track  disappearance  is  detected  when 
the  statistic  A*  =  L*_x  -  U*  exceeds  the  threshold  b  (the  lower  one).  When  the  decision  on 
target  disappearance  is  made,  the  statistic  U*  is  renewed  form  zero,  i.e.  the  detection  algorithm 
is  prepared  to  detect  another  target.  It  is  seen  that  the  algorithm  is  able  to  detect  even  very 
low  SNR  targets.  The  threshold  a  was  chosen  such  that  the  false  alarm  rate  1  false  alarm  per 
minute  (i.e.,  per  60  frames)  was  guaranteed.  The  algorithm  never  raised  more  than  1  false  alarm 
and  always  detected  targets  with  this  choice  in  the  analyzed  images.  To  be  precise,  there  are 
two  targets  in  the  pictures:  the  first  target  appears  at  time  n  =  1  and  disappears  at  n  =  28 
while  the  second  one  appears  at  time  n  =  39  and  does  not  disappear.  The  algorithm  detects  the 
first  target  with  the  delay  about  20  seconds  (20  frames)  for  SNR  =  — 6.6dB  and  4  —  6  frames 
for  SNR  =  — 0.25dB.  Then  the  fact  of  target’s  disappearance  is  detected  with  very  small  delay. 
Finally  the  second  target  is  detected  with  the  delay  about  5-6  frames  for  SNR  =  -6.6dB  and 
2  —  3  frames  for  SNR  =  — 0.25dB. 
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Adaptive  CUSUM  statistic,  no  target  Adaptive  CUSUM  statistic,  SNR=-0.25dB 


Figure  21.  Plots  of  the  adaptive  CUSUM.  Left  -  no  target.  Right  -  1st  target 
appears  at  n  =  1,  disappears  at  n  =  28;  2nd  appears  at  n  =  39,  doesn’t  disappear 


Adaptive  CUSUM  statistic,  no  target  Adaptive  CUSUM  statistic,  SNR=-3.62dB 


Figure  22.  Plots  of  the  adaptive  CUSUM.  Left  -  no  target.  Right  -  1st  target 
appears  at  n  =  1,  disappears  at  n  =  28;  2nd  appears  at  n  =  39,  doesn’t  disappear 


Adaptive  CUSUM  statistic,  no  target 


Adaptive  CUSUM  statistic,  SNR=-6.60dB 


Figure  23.  Plots  of  the  adaptive  CUSUM.  Left  -  no  target.  Right  -  1st  target 
appears  at  n  =  1,  disappears  at  n  —  28;  2nd  appears  at  n  =  39,  doesn’t  disappear 
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Algorithm  3,  SNR=-0.25dB 


Figure  24.  Detection  of  track  appearance  and  disappearance  by  Algorithm  3. 
First  target  appears  at  n  =  1  and  disappears  at  n  =  28;  second  target  appears  at 
n  =  39  and  doesn’t  disappear 
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8.  Conclusion.  Future  Research 

We  believe  that  real  time  optimal  nonlinear  filtering  is  an  emerging  information  technology  of 
fundamental  importance.  The  ONF  is  an  advance  over  conventional  EKF  techniques  and  multi¬ 
dimensional  (spatio-temporal)  matched  filters  potentially  as  profound  as  that  of  the  EKF  over 
a-j3-<y  filters.  Nonlinear  filtering  techniques  offer  the  promise  of  greatly  enhanced  performance 
in  a  number  of  DoD  related  applications:  airborne,  surface,  and  subsurface. 

We  have  developed  an  advanced  ONF-based  spatio-temporal  tracking  filter  and  adaptive  track 
appearance/disapperance  detection  algorithms  that  may  be  fully  incorporated  into  a  complete 
end-to-end  IRST  signal/track  processing  suite.  The  developed  technology  will  be  of  portable 
nature  and  can  be  incorporated  later  in  various  multisensor  ATR  systems  utilizing  EO/IR,  SAR, 
SAS,  LIDAR,  etc. 

In  the  future  we  plan  to  continue  this  work  in  the  following  key  directions. 

1.  Further  develop  and  tune  spatio-temporal  nonlinear  semiparametric  algorithms  for  clutter 
removal  and  jitter  compensation. 

2.  Develop  banks  of  ONF  filters  for  TbD  of  multiple  targets. 

3.  Improve  and  tune  track  detection  and  identification  algorithms. 

4.  Test  and  validate  the  developed  signal/data/track  processing  system  for  real  IR  data  (jointly 
with  Space  &  Naval  Warfare  Systems  Center,  San  Diego,  CA). 

5.  Incorporate  the  developed  system  into  a  complete  end-to-end  IRST  signal/track  processing 
suite  (jointly  with  Space  &  Naval  Warfare  Systems  Center,  San  Diego,  CA). 
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